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ABSTRACT 


Earth-Crossing  Asteroids  (EGAs)  are  those  asteroids  whose  orbit  cross-section 
can  intersect  the  capture  cross  section  of  the  Earth  as  a  result  of  secular  gravitational 
perturbations.  This  thesis  provides  a  framework  for  understanding  the  origin,  nature,  and 
types  of  EGAs.  The  change  in  velocity  requirements  to  achieve  a  two  Earth  radii 
deflection  for  long  and  short-term  warning  scenarios  are  developed.  Next,  a  method  of 
developing  hypothetical  Earth  colliding  asteroid  orbits  is  presented.  These  hypothetical 
orbits  are  used  in  two  ways:  (1)  to  evaluate  the  ability  of  Dance  of  the  Planets,  a  solar 
system  simulation  model  developed  by  Applied  Research  and  Consulting,  Inc.,  to 
accurately  propagate  orbits  of  imported  asteroid  orbits,  and  (2)  to  analyze  the  sensitivity 
of  deflection  distance  to  variation  in  deflection  angle  and  orbital  parameters  of  a  given 
orbit.  Inaccuracies  during  importation  of  data  precluded  the  use  of  Dance  of  the  Planets 
for  the  purpose  of  sensitivity  analysis.  The  program  does  provide  an  excellent  tool  for 
visualization  of  EGA  scenarios.  Consequently,  a  simpler  orbital  model  was  developed  to 
provide  a  Earth  miss  distance  sensitivity  analysis.  With  one  asteroid  orbital  period 
warning  the  minimum  change  in  velocity  to  deflect  an  asteroid  two  Earth  radii  is 
approximately  0. 135  m/s  and  the  optimal  deflection  is  along  the  flight  path.  Maximum 
deflection  occurs  when  the  deflection  is  applied  at  perihelion.  The  miss  distance  decreases 
markedly  with  increase  in  true  anomaly  until  it  is  a  minimum  at  aphelion. 
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I.  INTRODUCTION 


Impacts  by  Earth-crossing  Asteroids  (EGAs)  and  comets,  collectively  known  as 

Near-Earth  Objects  (NEOs),  pose  a  significant  and  unique  challenge  to  the  scientific 

community.  The  study  of  these  space-borne  bodies  has  spanned  a  broad  range  of 

disciplines:  mechanical,  electrical,  aeronautical,  and  astronautical  engineers, 

astronomers,  nuclear  physicists,  to  name  but  a  few.  Various  camps  have,  at  different 

times,  proposed  new  search  strategies  and  detection  methods  to  ensure  a  proper  tally  of 

potential  colliders,  forwarded  techniques,  both  nuclear  and  non-nuclear,  to  mitigate  the 

disaster  a  colliding  body  may  impart,  and  designed  missions  to  intercept,  rendezvous,  and 

study  a  few  of  the  closer  bodies.  The  recent  impact  of  comet  Shoemaker-Levy  9  with 

Jupiter  only  served  to  stimulate  additional  interest  in  the  subject  of  NEOs. 

Recognizing  the  potential  seriousness  of  such  events,  the  United  States  Congress 

in  1992  mandated  that  the  National  Aeronautics  and  Space  Administration  (NASA) 

conduct  two  workshops  to  study  two  major  research  areas  of  NEOs:  Detection  and 

Mitigation.  The  United  States  House  of  Representatives,  in  NASA  Multiyear 

Authorization  Act  of  1990  said,  in  part: 

The  committee  believes  that  it  is  imperative  that  the  detection  rate  of 
Earth-orbit-crossing  asteroids  must  be  increased  substantially,  and  that  the 
means  to  destroy  or  alter  the  orbits  of  asteroids  when  they  threaten 
collision  should  be  defined  and  agreed  upon  internationally. 

The  chances  of  the  Earth  being  struck  by  a  large  asteroid  are  extremely 
small,  but  since  the  consequences  of  such  a  collision  are  extremely  large, 
the  Committee  believes  it  is  only  prudent  to  assess  the  nature  of  the  threat 
and  prepare  to  deal  with  it.  We  have  the  technology  to  detect  such 
asteroids  and  to  prevent  their  collision. 

The  Committee  therefore  directs  that  NASA  undertake  two  workshop 
studies.  The  first  would  define  a  program  for  dramatically  increasing  the 
detection  rate  of  Earth-orbit-crossing  asteroids;  this  study  would  address 
the  costs,  schedule  technology,  and  equipment  required  for  precise 
definition  of  the  orbits  of  such  bodies.  The  second  study  would  define 
systems  and  technologies  to  alter  the  orbits  of  such  asteroids  or  to  destroy 
them  if  they  should  pose  a  danger  to  life  on  Earth.  The  Committee 
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recommends  international  participation  in  these  studies  and  suggest  that 

they  be  conducted  within  a  year  of  the  passage  of  this  legislation.  [Ref  1] 

As  a  result  two  conferences  were  conducted.  The  first,  the  NASA  International 
Near-Earth-Object  Detection  Workshop,  was  conducted  in  three  sessions  jfrom  June 
through  November  1991 .  Their  work  concentrated  on  improving  the  rate  at  which  EGAs 
are  discovered;  the  results  are  documented  in  reference  1.  The  second  conference.  The 
Near-Earth-Object  Interception  Workshop,  was  held  in  January  1992  and  hosted  by  Los 
Alamos  National  Laboratory.  It  concentrated  on  the  issues  surround  the  mitigation  of 
EGAs  and  associated  technologies.  Even  more  recently,  an  EGA  conference  was  held  in 
March  1995  at  Lawrence  Livermore  National  Laboratories.  Gonference  proceedings  have 
not  yet  been  published. 

These  conference  reports  address  the  majority  of  the  issues  surrounding  NEOs; 
however,  there  is  much  work  remaining.  Advances  in  detection  technology  such  as 
improved  sensors,  computer  search  algorithms,  d  space  radar  support  would 
significantly  speed  the  rate  at  which  NEOs  are  identified.  Space  missions  to  candidate 
asteroids  can  provide  essential  information  regarding  the  structure  and  composition  of 
asteroid  bodies.  Optimization  of  intercept  trajectories  and  continued  improvements  to  the 
space  launch  vehicles  will  all  be  directly  applicable  to  the  mitigation  of  these  potentially 
dangerous  objects. 

At  the  onset  of  this  project  there  were  three  main  goals:  (1)  Gonduct  a  thorough 
and  exhaustive  sun^ey  of  the  current  literature  in  the  area  of  EGAs,  (2)  Analyze  and 
select  fi-om  commercially  available  software  applications  the  candidate  most  likely  to  aid 
in  the  visualization  of  asteroid  mitigation  through  deflection,  and  (3)  Develop  mitigation 
scenarios  and  determine  the  sensitivity  of  the  circular  solution  to  elliptical  orbits. 

Research  under  the  Space  Warfare  Research  Program,  sponsored  by  the  Air 
Force’s  Space  Warfare  Genter,  was  divided  into  four  sections.  Two  major  components  of 
NEOs,  EGA’s  and  Earth-crossing  Gomets  (EGGs),  are  described  and  categorized. 
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Deflection  options  are  presented  and  the  maximum  size  asteroid  which  could  be  deflected 
is  estimated.  The  most  effective  deflection  angle  is  established  and  the  required  change 
in  velocity  to  perturb  an  object  by  two  Earth  radii  is  determined.  Methodology  for 
developing  hypothetical  asteroid  orbits  is  also  presented.  A  commercially  available 
software  application,  Dance  of  the  Planets,  is  evaluated  for  its  usefulness  in  verifying  the 
deflection  of  assailant  objects  and  testing  the  mathematical  solutions  for  optimal 
deflection.  Finally,  an  analysis  of  sensitivity  of  miss  distance  to  variations  in  orbit  shape 
and  deflection  direction  is  conducted.  The  issues  compiled  in  this  thesis  address  only  a 
small  portion  of  the  asteroid  mitigation  problem,  and  an  even  smaller  portion  of  the 
broader  subject  of  EGAs. 
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11.  NEAR  EARTH  OBJECTS 


A.  INTRODUCTION 

There  are  two  broad  categories  of  objects  whose  orbits  bring  them  close  to  that  of 
the  Earth:  asteroids  and  comets.  Astronomers  classify  the  objects  into  one  of  the  two 
categories  based  upon  their  telescopic  appearance.  Any  object  which  appears  to  be  star- 
like  is  called  an  asteroid.  If  the  object  has  a  visible  atmosphere  or  a  tail  then  it  is  a  comet. 
[Ref  1]  This  difference  is,  in  part,  due  to  the  composition  of  the  of  the  object.  Asteroids 
have  no  atmosphere  and  may  have  physical  and  compositional  properties  ranging  from 
loosely  aggregated  cometary  ices  to  solid  metal.  [Ref  2]  Comet  nuclei  are  composed  of 
a  complex  mixture  of  volatile  ices,  water,  and  hydrocarbon  and  silicate  grains.  As  a 
comet  is  heated  by  the  Sun  as  it  approaches  perihelion,  there  is  a  noticeable  out-gassing  of 
evaporative  material.  [Ref  3]  Further  discussions  on  each  broad  classification  of  Near- 
Earth  Objects  (NEOs)  are  provided  below. 

B.  ASTEROIDS 

When  viewing  the  solar  system,  shown  in  Figure  1  Figure  1,  the  empty  gap 
between  Mars  and  Jupiter  is  readily  apparent.  Kepler  contemplated  a  “missing  planet”  as 
he  studied  the  solar  system,  as  described  by  the  measurements  of  Tycho  Brahe.  [Ref.  4] 
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Figure  1.  Depiction  of  Solar  System  Through  Jupiter.  From  Ref.  4. 

Bode’s  law,  named  after  Johann  Bode’s  effort  to  stir  interest  in  an  unknown  planet 
in  1772,  uses  a  simple  relationship  to  generate  the  mean  distances  in  Astronomical  Units 
(AU)  of  the  principal  planets.  The  relationship  is  generated  by  writing  down  the  series  0, 
3, 6, 12, ...,  add  4  to  each  number  and  dividing  by  10.  As  shown  in  Table  1,  Bode’s  law 
generates  fairly  accurate  locations  for  all  the  planets  except  Neptune  and  Pluto. 


BODE’S  LAW 

ACTUAL 

PLANET 

DISTANCE 

DISTANCE 

Mercury 

0.4 

0.39 

Venus 

0.7 

0.72 

Earth 

1.0 

1.00 

Mars 

1.6 

1.52 

Asteroids  (average) 

2.8 

2.65 

Jupiter 

5.2 

5.20 

Saturn 

10.0 

9.52 

Uranus 

19.6 

19.28 

Neptune 

38.8 

30.17 

Pluto 

77.2 

39.76 

Table  1.  Bode’s  Law  vs.  Mean  Planetary  Distance 
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This  fits  the  actual  positions  of  the  planets  through  Saturn  very  well,  except  for  the 
position  at  2.8  AU.  This  position,  between  the  orbit  of  Mars  and  Jupiter  remains  empty, 
except  for  the  asteroid  belt.  [Ref.  4] 

1.  Early  History 

In  view  of  Bode’s  law,  a  group  of  24  astronomers  formed  a  society  in  Europe  in 
1 800  to  solve  the  problem  of  a  planet,  surmised  to  be  missing,  between  Mars  and  Jupiter. 
Each  astronomer  was  given  a  region  of  the  zodiac.  [Ref  4]. 

Giuseppe  Piazzi,  in  Palermo,  Italy,  was  already  engaged  in  a  star  charting  project. 
He  located  a  dim,  imcharted  ‘star’  that  shifted  position  fi'om  night  to  night.  This  was 
Ceres,  discovered  on  January  01,  1801.  He  observed  it  for  41  days  before  it  was  lost  due 
to  bad  weather  and  the  illness  of  Piazzi.  [Ref.  4]. 

Karl  Friedrich  Gauss  became  involved  with  the  intricate  mathematical  problem  of 
calculating  an  adequate  recovery  orbit  for  Ceres.  With  his  help,  Ceres  was  sighted  on 
December  07, 1801.  Gauss’  ingenuity  was  of  great  significance  to  the  field  of  celestial 
mechanics.  The  basic  elements  of  Ceres  were  eccentricity  of  0.08,  inclination  of  1 1 and 
semi-major  axis  of  2.77  AU.  [Ref.  4] . 

Heinrich  Wilhelm  Gibers  found  a  second  unknown  ‘planet’,  Pallas,  in  March 
1802,  while  helping  Gauss  observe  Ceres.  Gauss  next  calculated  an  orbit  for  Pallas. 

Many  scientists  contributed  in  the  calculation  of  the  orbital  elements,  and  in  the 
development  of  perturbation  theories  to  accoxmt  for  the  influence  of  Jupiter.  This  was 
one  of  the  first  documented  efforts  at  the  development  of  planetary  theory  from  which 
highly  accurate  planet  ephemerides  are  possible.  [Ref  4]. 

Discovery  of  additional  asteroids  progressed  rapidly  over  the  course  of  the  next 
century,  particularly  with  the  implementation  of  photography  in  the  search  process  in 
1891.  By  1900, 463  asteroids  had  been  discovered,  and  by  1950, 1568.  By  1993, 
approximately  5500  numbered  asteroids  had  been  documented.  Approximately  300 
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numbered  asteroids  are  added  annually.  Numbered  asteroids  have  orbits  confirmed  by 
observations  at  two  or  more  oppositions  (when  the  asteroid-Earth-Sun  angle  equals  180°). 
Many  others  asteroids  have  only  provisional  designations.  [Ref.  4]. 

2.  Origin  and  Nature  of  Asteroids 

There  are  many  plausible  explanations  for  the  existence  and  formation  of 
asteroids.  One  scenario  suggests  the  asteroids  formed  from  planetesimals  which  were 
never  able  to  accrete  into  planet  sized  bodies  due  to  the  disturbing  influences  of  a  proto- 
Jupiter.  There  are  many  inconsistencies  with  this  theory.  First,  from  the  planetesimal 
model,  there  should  be  between  one  and  two  Earth  masses  of  planetesimals  in  the  region 
of  the  asteroid  belt.  All  of  the  asteroids  together,  however,  are  no  more  than  0.08%  of  the 
mass  of  the  Earth.  Some  mechanism  is  responsible  for  this  loss  of  mass.  Secondly, 
individual  asteroids  have  high  relative  velocities.  They  do  not  fit  the  orderly  picture  of 
planetesimals  analogous  to  Saturn’s  rings:  bodies  lying  in  nearly  the  same  plane  with  low 
relative  velocities  between  neighboring  bodies.  Again,  some  mechanism  is  responsible 
for  the  high  relative  velocities  which  produce  collision  fragmentation  and  not  collision- 
accretion.  Finally,  many  stony  and  metallic  meteorites,  derived  from  asteroids,  have  clear 
indications  of  melting,  elemental  differentiation,  and  slow  cooling,  similar  to  that  foxmd 
in  the  planets.  The  process  of  this  development  is  not  fully  understood.  [Ref  4]. 

There  has  been  much  debate  as  to  the  origin  of  near-Earth  asteroids.  Because 
these  are  planet  crossing  asteroids,  they  have  a  dynamic  lifetime  (the  average  time  before 
a  planet  crossing  asteroid  impacts  one  of  the  inner  planets)  of  approximately  3  x  lO’  yr. 

As  such,  there  must  be  some  mechanism  responsible  for  their  replenishment.  There  are 
two  main  theories,  divided  largely  between  observers  and  dynamicists,  which  attempt  to 
explain  the  source  of  these  objects.  Observers  generally  believe  that  the  near-Earth 
asteroids  were  derived  from  the  main  asteroid  belt  due  to  spectropically  similar  objects  in 
the  main  belt.  Dynamicists  believe  they  may  be  largely  of  cometary  origin.  They  do  not 
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believe  that  there  is  sufficient  dynamical  interaction  in  the  main  asteroid  belt  to  resupply 
the  near-Earth  orbit  asteroid  population.  [Ref.  3] 

Asteroids  are  numbered  in  the  order  they  are  discovered  with  orbits  verified  with  a 
minimum  of  three  observations.  The  majority  also  have  names  assigned.  The  first 
asteroid  discovered,  1  Ceres,  is  approximately  932  km  in  diameter,  and  represents  about 
30%  of  the  total  mass  of  asteroids.  4  Vesta  and  2  Pallas,  approximately  528  km  and  523 
km,  respectively,  are  next  in  size.  There  are  approximately  30  other  main  belt  asteroids 
with  diameters  larger  than  200  km.  For  objects  with  the  same  density,  the  mass  varies  as 
the  third  power  of  their  linear  dimensions.  From  this,  it  is  easy  to  ascertain  that  most  of 
the  total  mass  of  the  asteroids  is  found  in  the  few  largest.  Ceres  is  just  below  the  limit  for 
the  class-3  satellite  size  category  of  1000-1600  km,  of  which  Saturn  and  Uranus  each 
have  four.  [Ref.  4]. 

3.  Asteroid  Close  Approaches 

Although  there  are  a  large  number  of  asteroids,  their  number  is  limited  by  our 
ability  to  detect  them.  This  remains  one  of  the  biggest  challenges  in  the  field.  Most 
newly  discovered  asteroids  are  dimmer  and  usually  smaller  than  previously  discovered 
asteroids.  By  far  the  majority  are  firagments  of  ancient  parent  bodies,  unlike  Ceres  which 
has  retained  most  of  its  original  shape  and  mass.  [Ref  4]. 

Near  Earth  Asteroids  (NEA)  are  of  interest  because  they  lie  well  inside  the  main 
asteroid  belt,  and  they  do  not  have  long  term  stability  because  they  will  eventually  either 
1)  collide  with  one  of  the  inner  planets,  or  2)  be  ejected  fi-om  this  region  by  a  near¬ 
collision  encounter.  Earth  impact  rate  for  objects  of  potentially  cataclysmic  size  is 
estimated  at  3-4  per  million  years.  No  such  collisions  have  been  yet  been  predicted, 
however.  The  “spacewatch”  telescope  on  Kitt  Peak  finds  several  close  approaches  per 
month.  [Ref.  4]. 
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Many  EGAs  have  come  relatively  close  to  the  Earth.  Hermes,  in  1937,  passed 
within  800,000  km  (twice  the  distance  to  the  Moon)  of  the  Earth.  It’s  orbit  was  lost  and 
it’s  location  is  no  longer  known.  On  January  18,  1991,  199  IB  A  passed  within 
approximately  170,000  km  (0.001 1  AU)  of  the  Earth.  The  newly  discovered  near-Earth 
asteroid  1994  XMl  made  the  closest  approach  to  the  Earth  of  any  object  discovered 
outside  the  earth's  atmosphere — some  105,000  km  on  the  morning  of  Dec.  9  over  Russia. 
The  diameter  was  estimated  to  be  approximately  6-13  meters. 

There  are  many  more  asteroids  with  the  potential  for  inner  planet  encounters. 
Among  them  are  433  Eros,  887  Alinda,  1036  Ganymed,  1221  Amor,  1566  Icarus,  1580 
Betulia,  1620  Geographos,  and  1627  Ivar.  There  are  many  other  inner  asteroids  of  note. 
3200  Phaethon  has  the  greatest  asteroid  eccentricity,  and  the  smallest  perihehon,  0. 135 
AU,  1951  Lick  is  located  entirely  between  the  orbits  of  Earth  and  Mars;  it  doesn’t  cross 
any  orbit.  2062  Aten  is  entirely  inside  Earth’s  orbit.  2102  Tantalus  has  the  highest 
inclination  of  any  numbered  asteroid  (64°).  1973na  has  an  even  greater  inclination  (68°). 

4.  Earth-Crossing  Asteroids 

There  are  many  asteroids  well  within  the  orbits  of  the  inner  planets.  A  current  list 
of  all  EGAs  through  1991  in  contained  in  reference  5.  These  asteroids  are  of  interest 
because  they  lie  well  inside  the  main  asteroid  belt.  Some  of  these  come  quite  close  to  the 
orbit  of  Earth.  An  Earth-crossing  Asteroid  (EGA)  has  been  defined  as  a  minor  planet 
whose  orbit  can  intersect  the  capture  cross-section  of  the  Earth  as  a  result  of  secular 
gravitational  perturbations.  [Ref.  5]  EGAs  have  secular  periods  on  the  order  of  tens  of 
thousands  of  years. 

a.  Secular  Variations 

It  has  been  shown  that  an  asteroid  orbit  which  overlaps  the  orbit  of  a  planet 
may,  as  a  result  of  the  advance  of  the  line  of  apsides,  intersect  the  orbit  of  the  planet.  If 
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the  orbit  of  the  planet  is  of  low  eccentricity  and  is  overlapped  by  an  asteroid  orbit  which  is 
highly  eccentric,  the  two  orbits  will  be  linked,  like  links  of  a  chain,  when  the  argument  of 
perihehon  is  0  or  tt  and  unlinked  when  the  unlinked  at  %I2  or  37r/2.  In  one  complete 
rotation  of  the  line  of  apsides,  the  orbits  transition  from  linked  to  unlinked  states  a  total  of 
four  times,  providing  four  possible  intersections  of  the  two  orbits.  These  intersections  of 
crossings  may  be  found  by  solving  the  polar  equation  of  the  ellipse,  and  occur  at 


Q  =  cos  '  + 


g(l-e^) 


1 


0) 


where  ®  is  the  argument  of  perihelion,  a  is  the  semi-major  axis  of  the  asteroid  orbit,  e  is 
the  eccentricity  of  asteroid  orbit  at  the  time  of  intersection,  and  p  is  the  radius  to  the 
planet’s  orbit  along  the  line  of  nodes  at  the  time  of  intersection.  [Ref  6]  Because  most 
planet’s  orbits  are  not  circular  and  there  are  secular  variations  'm  e,  p  will  have  four 
different  values  for  each  rotation  of  the  line  of  apsides  which  can  be  found  by 
simultaneously  solving  Eqn  (1)  and  the  polar  equation  for  the  elliptical  orbit  of  the  planet. 
This  type  of  orbit  is  referred  to  as  a  quadruple  crasser,  and  is  shown  in  Figure  2  using 
asteroid  2062  Aten  as  an  example.  The  figure  shows  one  cycle  of  advance  of  ©  starting  at 
©  =  0.  The  solid  line  shows  the  ascending  node  and  the  dashed  line  shows  the  descending 
node. 


Figure  2.  Secular  Variation  of  Radius  to  the  Node  of  the  Orbit  of 2062 
Aten.  From  Ref  6. 
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Secular  variations  in  e  may  be  sufficiently  high  to  change  the  depth  of 
orbital  penetration  of  the  assailant  body.  In  this  case  the  orbits  may  transition  from  a  non- 
overlapped  condition  at  co  =  0,  tt  to  a  relatively  deep  overlap  at  co7c/2, 3nl2.  The  first 
crossing  occurs  due  to  the  secular  increase  in  e  as  co  begins  its  rotation  from  0  to  %I2,  thus 
linking  the  orbits.  The  second  crossing  occurs  as  co  increases  sufficiently  to  unlink  the 
orbits.  Thus,  there  are  two  crossings  during  each  V2  advance  of  the  line  of  apsides.  This 
type  of  asteroid  is  known  as  an  octup:  crasser.  An  example  is  shown  in  using  asteroid 
1580  Betulia,  the  first  asteroid  where  this  type  of  behavior  was  recognized.  [Ref.  6] 
Approximately  1.5  cycles  of  advance  of  co  are  shown,  viith  the  solid  line  representing  the 
ascending  node  and  the  dashed  line  representing  the  descending  node. 


Figure  3.  Secular  Variation  of  Radius  to  the  Node  of  the  Orbit  of 
1974MA.  From  Ref.  6. 

If  the  inclination  of  the  asteroid  with  respect  to  the  ecliptic  plane  is 
sufficiently  high,  the  secular  perturbations  are  not  strong  enough  to  influence  a 
continuous  advance  in  co.  In  this  case,  co  librates  around  tcI2  or  3%I2.  When  this  occurs,  a 
large  oscillation  of  e  can  occur  during  the  libration  cycle  leading  to  orbit  intersection  four 
times  during  each  libration  cycle  of  co.  These  types  of  asteroid  are  called  quadruple 
crossing  co  librators.  An  example  is  shown  in  Figure  4  using  asteroid  1973  NA. 
Approximately  five  libration  cycles  of  omega  are  shown.  [Ref  6}. 
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Figure  4.  Secular  Variation  of  Radius  to  the  Node  of  the 
Orbit  of  1973  NA.  From  Ref.  6. 


The  fourth  type  of  EGA  are  known  as  supercrossers  and  result  from 
asteroids  that  librate  about  the  3:1  commensurability  with  Jupiter.  The  resonant  effects  of 
these  perturbations  cause  a  relatively  high  frequency  oscillation  in  a  and  e  which  lead  to 
relatively  high  frequency  oscillations  (4  cycles  in  approximately  1400  years)  between 
periods  of  overlap  and  nonoverlap  of  the  two  orbits.  For  example,  asteroid  1915 
Quetzalcoatl,  shown  in  Figure  5,  had  nine  crossings  of  the  Earth’s  orbit  over  a  period  of 
approximately  1400  years.  Approximately  four  cycles  of  libration  of  the  mean  motion  of 
the  asteroid  are  shown.  [Ref  6]. 


Figure  5 .  Secular  V ariation  of  Radius  to  the  Ascending 
Node  of  1915  Quetzalcoatl.  From  Ref.  6. 

b.  Classes  of  Earth-Crossing  Asteroids 

In  addition  to  classifying  asteroids  by  the  type  of  secular  perturbations  they 
experience,  Earth-crossing  asteroids  are  further  divided  into  three  groups  on  the  basis  of 
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their  present  osculating  orbital  elements:  (1)  Aten  asteroids,  (2)  Apollo  asteroids,  and 
(3)  Earth-crossing  Amor  asteroids. 

(1)  Aten.  The  orbits  of  Aten  asteroids  have  a  semi-major  axis  less 

than  1  AU.  Their  orbits  overlap  those  of  the  Earth  near  their  aphelia,  as  shown  in  Figure 

6. 


View  Looking  Down  +Z 
.  Y 


Asteroids  Orbit 


2340  Hathor 
a  =  0.844 
e  =  0.450 


Figure  6.  An  Aten  Type  Asteroid. 

This  class  includes  all  asteroids  with  a  <  1.0  AU  and  with  the  aphelion  distance  of  the 
asteroid,  Q  >  0.983  AU.  The  asteroid’s  orbit  can  intersect  that  of  Earth  as  the  Earth’s 
perihelion  distance  is  0.983  AU.  The  first  Aten  was  discovered  relatively  recently,  in 
1976.  Fifteen  were  known  as  of  August  1993.  They  represent  approximately  10%  of  the 
180  known  EGAs.  The  first  three  discovered  exhibit  continuous,  or  nearly  continuous 
orbital  overlap  with  Earth  and  are  characterized  as  deep  quadruple  crossers.  [Ref.  6]. 

The  total  number  of  Aten  asteroids  to  visual  magnitude  has  been 
estimated  to  be  on  the  order  of  100.  In  the  same  way  that  some  Earth-crossing  Amors 
have  current  osculating  orbits  entirely  outside  the  orbit  of  the  Earth,  some  Atens  may 
have  ort  Its  entirely  inside  the  orbit  of  the  Earth.  Shoemaker  et.  al.  numbers  these  at  a  few 
tens  of  objects  out  to  visual  magnitude  18.  [Ref  6] 
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(2)  Apollo.  The  orbits  of  Apollo  asteroids  are  those  asteroids  with 
a>  1.0  AU,  and  perihelion  distance,  ^<1.017  AU,  where  1.017  is  the  aphelion  distance  of 
the  Earth.  An  example  is  presented  in  Figure  7. 

View  Looking  Down  *Z  axis 


Figure  7.  An  Apollo  Type  Asteroid. 


125  Apollos  have  been  discovered,  which  represents  approximately  two-thirds  of  all 
EGAs.  All  but  about  5%  of  all  Apollo  asteroids  are  EGAs.  [Refs.  6  and  5]  Of  those 
known  in  1979,  approximately  60%  were  quadruple  crossers,  20%  were  octuple  crossers, 
5%  were  quadruple  crossers  part  of  the  time  at  octuple  crossers  part  of  the  time,  and  15% 
were  quadruple  crossing  ©  librators.  [Ref.  6]  In  1932  the  first  Earth-crossing  asteroid, 

1 862  Apollo,  was  found  by  K.  Reinmuth  at  Heidelberg.  [Ref.  8] 

(3)  Amor.  Although  Amor  asteroids  have  perhaps  the  simplest 
definition  of  all,  the  fact  that  some  are  in  Earth-crossing  orbits  can  be  confusing.  Their 
orbits  are  defined  strictly  upon  the  orbital  perihelion  distance,  1.017  <  ^  <  1.3  AU.  They 
have  perihelions  close,  but  a  little  outside  of  the  Earth’s  orbit.  The  first  discovery  of  an 
asteroid  of  this  class  was  433  Eros  in  1 898.  [Ref.  4]  However,  the  traditional  name  for 
Atens  comes  firom  the  asteroid  of  the  same  name  discovered  by  E.  Delporte  at  Uccle,  in 
1932.  [Ref  8]  The  47  known  Amor  asteroids  comprise  approximately  25%  of  all  EGAs. 
[Ref.  5]  The  upper  limit  for  perihelion  distance  of  1.3  AU  was  chosen  because  is  it  is 
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near  a  minimum  in  the  radial  frequency  distribution  of  q  for  discovered  objects.  [Ref  6] 
As  discussed  previously,  some  Amors  can  become  Apollos  and  vice  versa,  due  to  secular 
perturbations.  The  classifications  are  based  on  the  category  at  the  time  of  discovery.  For 
example,  if  the  Amor  asteroid  1915  Quetzalcoatl  had  been  discovered  in  1942,  it  would 
have  been  classified  as  an  Apollo  asteroid.  [Ref  6]  An  example  of  an  Earth-crossing 
Amor  is  shown  Figure  8. 


View  Looking  Down  +Z  axis 


Figures.  An  Amor  Type  Asteroid. 

c.  Distant  Asteroids 

While  most  asteroids  are  contained  within  the  main  belt,  there  are  some 
that  lie  outside  that  area.  The  Trojans,  shown  in  Figure  1,  are  a  group  of  approximately 
100  asteroids  in  two  thinly  populated  lobes  centered  on  Jupiter’s  orbit  at  about  the 
Lagrangian  L4  and  L5  points.  They  are  in  stable  orbits  and  have  an  orbital  resonance  of 
1:1  (the  ratio  of  the  number  of  orbits  of  the  asteroid  to  that  of  a  third  body)  with  Jupiter. 
While  not  tightly  grouped  about  the  Lagrangian  points,  they  librate  about  these  nominal 
positions,  moving  closer  and  further  from  Jupiter,  but  not  far  from  its  orbit.  [Ref  4] 
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Opposite  Jupiter,  there  is  a  thinner  lobe  centered  on  the  Lagrangian  L3 
point.  These  are  known  as  the  Hilda  asteroids  and  have  a  3:2  resonance  with  Jupiter. 
Over  40  Hilda  type  asteroids  have  been  discovered.  Hilda  orbits  are  very  similar  to  some 
short-period  comet  orbits.  The  difference  being  that  the  Hilda  asteroids  are  in  stable 
orbits  while  the  short-period  comets  are  not.  [Ref.  4] 

Another  asteroid  of  note  is  279  Thule.  It  is  the  only  known  4:3  asteroid.  It 
is  unique  in  that  having  a  resonance  near  unity  places  it  near  Jupiter’s  orbit  where  it  is 
strongly  perturbed  every  three  Jovian  years,  yet  remains  in  a  stable  orbit.  944  Hidalgo  has 
a  comet-like  orbit  that  takes  it  out  nearly  to  Saturn.  1373  Cincinnati  experiences 
relatively  strong  perturbations  from  Jupiter  but  has  no  resonance  making  it  unique. 

2060  Chiron  was  the  first  trans-Satum  asteroid.  Discovered  in  1977,  it  has 
a  diameter  of  approximately  200  km.  Its  motion  has  been  simulated  by  researchers  back 
to  1664  BC,  when  it  appears  to  have  come  within  approximately  0.1  AU  of  Saturn. 
1992QB1  appears  to  be  of  similar  size  to  Chiron  but  its  orbit  extends  beyond  that  of 
Pluto.  Exotic  bodies  such  as  these  may  be  from  the  long-postulated  Kuiper  belt, 
discussed  in  the  comet  section. 

d.  Asteroid  Classifications 

Asteroids  are  of  four  different  classifications  based  upon  the  photometric 
characteristics  of  the  reflected  light,  the  amount  of  which  is  a  function  of  their  nature  and 
composition  of  their  surface  material.  This  classification  is  only  based  upon  spectral 
characteristics.  The  main  types  have  good  correlation  to  location  in  the  asteroid  belt. 

Most  asteroids  fall  into  one  of  the  following  four  categories.  [Ref  4]. 

S-type.  A  broad  distribution  of  asteroids  centered  at  about  2.3  AU 
from  the  Sun.  This  type  is  associated  with  stony-iron  meteorites,  although  some  scientists 
question  this  association.  These  asteroids  have  high  concentrations  of  nickel  and  iron, 
with  approximately  equal  amounts  of  metals  and  silicates.  These  are  presumed  to  be  the 
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exposed  fragments  from  the  inner  cores  of  larger  parent  bodies.  These  are  the  second 
most  prevalent  type  of  asteroid.  [Ref.  4]. 

C-type.  The  most  common  type  of  asteroid  with  a  relatively  sharp 
distribution  centered  at  about  3.1  AU  from  the  Sun.  Thought  to  be  composed  of  material 
similar  to  carbonaceous  chondritic  meteorites,  C-types  are  very  dark  with  an  albedo  of 
approximately  0.04.  They  have  various  grainy,  carbon  rich,  rock-like  mineralogies,  and 
are  a  complicated  group  to  examine  photospectfally.  Most  scientists  agree  Ceres  is  a  C- 
type  asteroid.  [Ref  4]. 

M-type.  Asteroids  with  moderate  reflectivity  and  a  reddish  spectra 
associated  with  metals,  particularly  Ni-Fe  are  known  as  M-type  asteroids.  Similar 
meteorites  are  the  nickel-irons  and  the  enstatite  chondrites,  which  have  Ni-Fe  embedded 
in  silicates.  There  are  not  many  M-type  asteroids.  The  distinction  between  M-type  and 
S-type  asteroids  is  not  clear  and  there  are  many  similarities  between  them.  [Ref.  4]. 

D-type.  D-type  asteroids  have  low  albedos  and  reddish  spectra  and 
are  primarily  the  Trojans,  located  near  the  orbit  of  Jupiter.  They  are  composed  mostly  of 
clays  with  magnetite  and  carbon-rich  materials.  No  meteorites  match  these  spectral 
characteristics  but  the  dark  side  of  Saturn’s  lapetus  is  a  good  match.  [Ref  4]. 

e.  Earth-Crossing  Asteroid  Population 

Visual  magnitude  is  defined  as  an  arbitrary  number,  measured  on  a 
logarithmic  scale,  used  to  indicate  the  brightness  of  an  object.  A  one  magnitude 
difference  is  the  fifth  root  of  100,  and  is  approximately  equal  to  a  factor  of  2.512.  The 
brighter  the  asteroid,  the  lower  the  numerical  value  of  magnitude.  Very  bright  objects 
have  negative  magnitudes. 

The  absolute  magnitude,  H,  of  an  asteroid  is  the  magnitude  the  asteroid 
would  have  if  it  were  1  AU  from  the  Sun,  and  viewed  from  the  Sim.  For  example,  if  an 
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asteroid  was  2  AU  from  the  Sun  and  viewed  from  1  AU  away,  it  would  be  four  times 
dimmer,  or  about  1.5  magnitude  less  bright. 

As  most  EGAs  are  discovered  through  direct  telescopic  observation,  the 
absolute  magnitude  of  the  asteroids  is  of  extreme  importance  to  an  observer.  In  general, 
the  larger  the  diameter  of  an  object,  the  lower  the  value  of  absolute  magnitude  (the  object 
is  brighter  than  a  similar  asteroid  of  smaller  diameter).  As  a  result,  the  large  asteroids  are 
often  discovered  more  quickly  than  the  smaller  objects. 

The  two  most  common  types  of  asteroids  are  C-type  and  S-type.  Table  2 
presents  the  current  estimates  of  EGA  discovery  completeness  based  on  absolute 
magnitude,  H.  Estimated  diameters,  based  on  the  albedo  of  the  two  most  prevalent  types 
of  asteroids,  are  also  presented. 


H 

Percent 

discovered 

Asteroid  type 

G-type  diameter 

13.2 

100% 

12  km 

6  km 

15.0 

35% 

6  km 

3  km 

16.0 

15% 

4  km 

2  km 

17.7 

7% 

2  km 

1  km 

Table  2.  Estimated  Discovery  Gompleteness  From  Ref.  5 


Numerous  studies  have  attempted  to  estimate  the  total  number  of  asteroids 
in  Earth-crossing  orbits.  These  estimates  are  based  on  a  power  law  relationship 

N  =  kD'’  (2) 

where  N  is  the  number  of  asteroids  larger  than  a  given  diameter,  D.  A:  is  a  constant  and  b 
is  the  power-law  exponent.  The  general  form  of  the  size  distribution  is  based  on 
observation  as  the  actual  detailed  distributions  remain  unknown. 

There  are  approximately  180  known  Earth-crossing  asteroids.  The 


characteristics  of  Apollo  asteroids  that  overlap  the  Earth’s  orbit  part  of  the  time  and  the 
characteristics  of  Earth-crossing  Amors  are  essentially  identical.  The  distinction  is 


derived  from  the  present  state  of  the  cycle  of  variation  of  their  perihelion  distance.  The 
majority  of  Earth-crossing  Amors  are  shallow  crossers,  while  the  majority  of  Apollos  are 
deep  crossers.  Because  the  orbits  of  most  Earth-crossing  Amors  only  overlap  the  orbit  of 
Earth  for  a  small  period  of  time,  they  retain  their  traditional  classification  as  Amor 
asteroids.  [Ref  6] 

Based  on  the  best  information  to  date,  an  estimate  of  the  number  of  Earth¬ 
crossing  asteroids  larger  than  a  given  diameter,  D,  are  shown  in  Figure  9.  For  this 
population  model,  changes  in  the  power  law  are  estimated  to  occur  at  diameters  of  10  m, 
70  m,  and  3.5  km. 


Figure  9.  Estimated  Number  of  Earth-Crossing  Asteroids.  From  Ref  5. 

C.  COMETS 

Comets  are  a  diffuse  bodies  of  gas  and  solid  particles  (such  as  CN,  C2,  NH3,  and 
OH),  which  orbit  the  Sun.  Their  orbits  are  highly  elliptical  or  even  parabolic  in  nature. 
Edmund  Halley  used  Newton’s  celestial  mechanics  to  show  that  comets  orbited  the  Sun. 
He  then  deduced  that  one  spectacular  comet  in  particular  was  returning  to  Earth  with  a 
period  of  76  years.  His  correct  prediction  of  the  comet’s  return  in  1758  made  it  comet 
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Halley.  Today,  as  many  as  lO"  to  lO’^  comets  remain,  orbiting  the  solar  system  in  the 
Oort  Cloud  and  Kuiper  Belt,  described  in  Section  2,  below. 

1.  Comet  Classifications 

Comets  which  have  an  elliptical  orbit  are  known  as  periodic  comets.  A  significant 
number  of  comets  also  have  parabolic  or  hyperbolic  orbits  and  will,  therefore,  make  only 
one  perihelion  passage.  The  periodic  comets  are  further  divided  into  short-period  and 
long-period  comets. 

a.  Short-Period  Comets 

Short  period  comets  are  defined  as  those  whose  orbits  lie  predominately 
within  the  realm  of  the  solar  system.  These  include  all  the  comets  which  have  been  seen 
more  than  once.  The  shortest  known  period  is  that  of  comet  Encke  at  3.3  years  while  the 
longest  extend  up  to  200  years.  Most  short-period  comets  have  inclinations  close  to  the 
ecliptic,  less  than  approximately  35°.  [Refs.  4  and  3]  About  95%  move  in  a  direct  sense. 
Comet  Halley  is  a  short  period  comet  in  a  retrograde  orbit.  More  recently,  short-period 
comets  have  been  further  divided  into  two  additional  classes:  (1)  Jupiter  and  (2)  Halley. 
Comets  in  the  Jupiter  family  have  periods  of  less  than  20  years  with  inclinations  close  to 
that  of  the  ecliptic  and  direct  orbits.  Halley  type  comets  have  period  of  20  <  P  <  200 
years.  In  general,  Halley  type  comets  have  higher  average  inclinations  than  Jupiter  type 
comets. 


b.  Long-Period  Comets 

Long  period  comets  are  those  with  periods  greater  than  200  years.  Long- 
period  comets  are  randomly  distributed  on  the  celestial  sphere.  Often  it  is  difficult  to 
distinguish  an  orbit  with  a  very  long  period  (on  the  order  of  10^  -  lO’  years)  from  a 
parabolic  or  hyperbolic  orbit.  Many  appear  to  have  aphelia  of 20,000-50,000  AU.  These 
types  of  comets  come  from  all  directions.  There  is  no  relationship  for  the  orientation  of 
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the  orbit  with  respect  to  the  ecliptic  or  for  direct  revolution.  Approximately  one-third  of 
all  long-period  comets  observed  are  on  weakly  hyperbolic  orbits. 

2.  The  Origin  and  Nature  of  Comets 

Many  scientists  believe  comets  to  be  planetesimals  formed  in  the  colder  reaches 
of  the  solar  system,  beyond  Saturn,  and  perhaps  beyond  Neptune.  Those  “cometesimals” 
in  the  planet  realm  of  the  solar  system  were  presumed  to  be  ejected  by  perturbations  from 
the  outer  planets.  These  comets  now  reside  in  the  Oort  Cloud,  at  the  outer  reaches  of  our 
Sun’s  gravitational  influence.  The  Oort  cloud,  named  after  Jan  Oort  in  1950,  is  a 
reservoir  of  comets  located  about  a  light-year  from  the  Sun. 

Comets  are  brought  into  the  Earth’s  realm  from  their  orbits  in  the  Oort  cloud  or 
the  Kuiper  belt  as  they  are  perturbed  by  Jupiter  or  passing  stars  every  few  million  years 
and  result  in  comet  insertion  back  into  the  solar  system.  [Ref.  9] 

At  the  same  time  Jan  Oort  was  formulating  his  hypothesis,  the  structure  of  comets 
was  proposed  by  Fred  Whipple.  His  “dirty  snowball”  model  hypothesizes  a  composition 
of  various  ices,  predominantly  H2O,  along  with  lesser  amounts  of  silicate  and  other 
mineral  dusts.  This  analysis  of  the  dirty  snowball/Oort  cloud  has  remained  the  generally 
accepted  model,  with  subsequent  modification  and  elaboration  due  to  follow-on  research. 
Relatively  recent  comet  space  probes,  particularly  to  comet  Halley  in  1986,  have  greatly 
increased  general  knowledge  about  the  comet. 

3.  The  Oort  Cloud 

The  factors  that  led  to  Oort’s  comet  cloud  model  were  1)  very  long-period  orbits 
that  were  commonly  found  to  have  aphelia  of  around  50,000  AU,  2)  the  large  reservoir 
that  would  be  required  in  order  to  supply  comets  over  the  last  4  billion  years  at  the  rate 
they  appear  in  the  inner  solar  system  (a  trillion  or  so),  and  3)  the  fact  that  very  long- 
period  comets  appear  from  all  directions.  Oort  proposed  that  the  planetary  system  was 
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siOTOunded  by  a  distant  spherical  cloud  of  comets  extending  approximately  one-half  the 
distance  to  the  nearest  stars.  The  comet  cloud  is  perturbed  by  passing  stars  which  cause 
the  injection  of  a  long-period  comet  into  our  solar  system.  The  source  of  the  Oort  cloud 
is  postulated  to  be  icy  planetesimals  ejected  by  the  proto-planets  in  the  outer  solar  system, 
particularly  Uranus  and  Neptune.  [Ref.  3] 

4.  The  Kuiper  Belt 

The  Kuiper  Belt  is  a  disk  thought  to  be  composed  icy  planetesimal  remnants 
beyond  Neptune,  and  was  first  proposed  by  Kuiper  in  1951.  Because  these  remnants  have 
such  a  high  orbital  periods  and  because  the  density  of  material  in  the  solar  nebula 
accretion  disk  decreases  beyond  Neptune,  this  material  was  unable  to  accrete  into  a 
planetary  body.  It  has  been  postulated  that  the  Kuiper  Belt  is  responsible  for  supplying 
the  low-inclination  Jupiter-class  comets.  [Ref.  3] 
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III.  DEFENDING  THE  EARTH  AGAINST  IMPACTS  FROM 
LARGE  COMETS  AND  ASTEROIDS. 


A.  INTRODUCTION 

There  are  numerous  illustrations  of  asteroid  and  comet  impacts  on  Earth: 
meteorite  craters,  astroblemns,  and  certain  circular  geologic  structures.  Many  hypotheses 
have  been  formulated  to  predict  the  consequences  of  a  collision  with  a  Near  Earth  Object 
(NEO)  with  the  Earth.  Damage  could  be  localized  to  a  small  geographical  area,  result  in 
regional  political  instabilities,  or  cause  devastating  climatic  change,  dwarfing  the  feared 
“nuclear  winter”  damage  of  a  global  thermonuclear  war.  Recent  research  into  the  effects 
of  a  collision  have  fairly  well  established  the  evidence  for  severe  climate  alterations.  The 
best  known  and  most  extensively  investigated  impact  occurred  on  the  Yucatan  Peninsula, 
at  the  boundary  of  the  Carboniferous  and  Paleogeneous  Ages,  65  million  years  ago,  and  is 
thought  to  have  caused  the  extinction  of  giant  reptiles  of  that  time.  The  geochemical  and 
paleontological  record  has  demonstrated  that  a  10-  to  15-km  diameter  NEO  impacted  the 
Earth  with  the  force  of  100  million  megatons  in  explosive  energy.  [Ref  1]  The 
frequency  of  such  impacts  is  small  but  not  insignificant. 

The  most  energetic  event  of  this  century  occurred  over  the  Tunguska  Valley  in 
1908.  It  was  surmised  that  this  event  resulted  from  the  break-up  of  a  chondritic  asteroid 
at  an  altitude  of  several  kilometers.  The  explosive  yield  was  originally  estimated  at  15- 
20  megatons  (Mt)  of  TNT,  although  recent  estimates  have  placed  that  number  as  high  as 
48  Mt,  which  leveled  approximately  2,000  square  kilometers  of  forest.  [Ref.  9]. 

There  are  numerous  ways  of  deflecting  an  Earth-crossing  asteroid  or  comet.  One 
involves  a  direct,  kinetic  energy  type  weapon,  which  imparts  a  change  in  momentum  to 
the  assailant  object.  A  second  way  of  mitigating  the  threat  from  an  EGA  or  comet  is 
through  the  use  of  nuclear  explosive  devices.  Issues  of  great  importance  surrounding  the 
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use  of  a  nuclear  device  on  a  potential  colliding  object  include;  1)  estimating  the 
radiation  energy  transfer  to  its  surface;  2)  maximizing  radiation  energy  transfer;  3) 
understanding  the  hydrodynamic  motion  and  mass  blow-off  from  the  surface,  and  the 
resultant  momentum  transfer;  and  4)  gaining  knowledge  of  the  various  ways  the  object 
may  be  disrupted  due  to  fracture  or  fragmentation.  [Ref  9],  Regardless  of  whether  one 
of  these,  or  some  other  method,  is  chosen,  the  several  basic  issues  must  be  thoroughly 
understood  before  an  attempt  is  undertaken  to  deflect  an  object. 

The  most  nkely  candidate  for  mitigating  the  danger  of  an  Earth  impact  is  a  rocket- 
delivered  nuclear  explosive.  Nuclear  explosive  devices  possess  the  highest  concentration 
of  energy,  and  can  be  manufactured  with  warheads  yielding  100  Mt  or  more.  A 
thermonuclear  charge  with  a  yield  of  1  Mt  has  a  mass  of  the  order  of  0.5  ton.  Of  interest 
is  the  fact  that  to  obtain  the  equivalent  energy  by  impact  of  a  body  of  the  same  mass 
requires  an  impact  velocity  of  approximately  4,000  km/s.  [Ref  9]. 

Simonenko  et  al.  [Ref  9]  presents  a  list  of  problems  that  must  be  considered 
before  the  development  of  a  defensive  system  could  be  undertaken. 

1 .  Assessment  of  the  number  of  dangerous  cosmic  bodies  and  their  probability  of 
their  collision  with  the  Earth. 

2.  Assessment  of  the  time  required  to  detect  and  identify  them  as  dangerous. 

3.  Understanding  of  the  object’s  motion  in  the  immediate  vicinity  of  the  Earth, 
penetration  of  the  atmosphere,  and  detailed  dynamics  during  collision  with  the  Earth, 
with  concomitant  assessment  of  the  local  and  global  consequences  of  the  collision. 

4.  Assessment  of  the  required  effect  on  the  astral  assailant,  whether  it  be 
deflection  or  fragmentation. 

5.  Consideration  of  the  needed  nuclear  devices,  delivery  means,  and  optimal 
regime  of  action. 

6.  Assessment  of  consequences  of  collision  with  fragments  of  an  object  that  had 
been  fractured  by  a  nuclear  explosive. 
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7.  Consideration  of  ecological  consequences  for  the  Earth  and  space  environment 
of  nuclear  explosions  in  space. 

The  most  important  question  from  a  technical  viewpoint  is  to  ascertain  if  it  is 
feasible  to  significantly  alter  the  trajectory  of  an  EC  A  or  comet  by  detonating  a  nuclear 
explosive  near  the  surface.  Equally  as  important  is  the  ability  to  calculate  the  effects  of 
such  a  weapon  on  the  extraterrestrial  body.  Many  characteristics  of  nuclear  explosions 
are  predictable  and  well  defined.  For  example,  detonation  of  a  100-kt  device  known  as 
“Sedan”  produces  a  crater  370  m  in  diameter  and  100  m  in  depth.  Significant  differences 
may  result  by  varying  the  nuclear  explosion  yield,  distance  to  the  target  object,  the  object 
dimensions,  material  composition  of  the  object,  and  the  design  of  the  nuclear  device. 
Depending  upon  the  interplay  of  these  variables,  the  resultant  explosion  can  result  in 
fragmentation,  crushing,  or  deviation  of  the  target  object  from  its  initial  trajectory.  [Ref. 

9]. 

Table  3  presents  the  yield  vs.  mass  for  nuclear  explosive  devices.  This  shows  that 
the  yield  from  current  nuclear  technology,  that  which  can  be  launched  with  existing 
rocket  technology,  is  capable  of  deflecting  small  objects.  The  nuclear  device  can  be 
scaled  upwards,  perhaps  an  order  of  magnitude  or  more,  while  preserving  special 
characteristics.  In  this  case,  modification  would  not  require  further  testing.  [Ref  9]  Of 
particular  importance  to  the  deflection  of  an  asteroid  are  the  specific  physical 
characteristics  of  the  body.  These  include  shape,  material  composition,  and  physical 
stability  of  the  object.  Maximum  effectiveness  for  deflection  may  require  several  nuclear 
devices  detonated  in  a  carefully  determined  sequence,  or  one  device  of  special 
configuration.  In  addition,  special  measures  can  be  taken  to  minimize  the  radiation 
contamination  of  the  asteroid  fragments,  particularly  important  if  it  is  determined  that  the 
fragments  from  such  an  explosion  will  continue  towards  impact  with  Earth. 
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Yield 

Mass 

1  Mt 

0.5  ton 

10  Mt 

3  to  4  ton 

100  Mt 

20  to  25  ton 

Table  3.  Yield  vs.  Mass  for  Nuclear  Devices  [Ref.  9] 


Finally,  an  area  of  study  which  needs  additional  study  is  that  of  the  explosion  of 
nuclear  devices  in  space.  Major  differences  include  absence  of  atmosphere, 
commensurability  of  the  object’s  dimensions  with  linear  scales  of  the  phenomenon, 
complexity  of  object  shape,  relatively  weak  gravity,  and  exotic  material  composition. 
[Ref  9] 

B.  DEFLECTION  CONSIDERATIONS 

There  are  two  main  areas  of  interest  in  the  mitigation  of  a  Near-Earth  Object 
(NEO):  1)  causing  the  threatening  object  to  break  up,  and  2)  altering  the  trajectory  of 
the  object.  There  are  several  issues  surrounding  these  areas  which  concern  the  scientific 
community  with  regard  to  deflecting  or  fragmenting  a  NEO  sufficiently  to  avert  collision 
with  Earth.  First,  is  there  a  limitation  to  the  size  of  asteroid  that  can  be  deflected,  given 
today’s  technology.  Secondly,  given  the  option,  what  is  the  most  effective  way  to  deflect 
these  astral  assailants  to  ensure,  with  a  reasonable  degree  of  certainty,  that  they  will  miss 
the  Earth. 

1.  Size  Limitations 

Simonenko  et  al.  [Ref.  9]  has  estimated  the  maximum  dimensions  of  an  object 
that  can  be  influenced  by  a  nuclear  explosion  as  follows.  Assume  a  nuclear  explosion 
occurs  on  the  surface  of  an  asteroid  having  mass  ma  and  radius  Ra.  Assume  the  mass  m 
of  the  material  ejected  by  the  explosion  is  small  compared  with  the  mass  of  the  object. 
The  mass  of  material  ejected  from  the  body  has  a  distribution  of  velocities,  but  to 
simplify  the  estimation,  it  is  assumed  that  all  the  ejected  material  has  the  same  velocity, 


28 


the  average  of  all  the  actual  velocities,  v„ .  This  varies  depending  upon  the  explosive 
yield  and  material  properties  of  the  object. 

Assuming  m«ma ,  conservation  of  momentum  shows 

mv„  =  (3) 

where  v«  is  the  velocity  of  the  ejected  material  an  infinity  and  Av^  is  the  velocity  increase 
of  the  object  as  a  result  of  the  explosion.  In  similar  fashion,  conservation  of  energy 
shows 

R,  2  ~  2  ^  2 


where  G  is  the  universal  gravitational  constant.  Combining  Eqns.  (3)  and  (4)  gives 


Vco  =  v„ 


2Gm„ 

.  ^  — o 


(5) 


vIr. 


Examination  of  Eqn.  (5)  shows  that  if  an  asteroid  is  of  sufficient  mass,  and  given 
the  above  limitations  on  ejected  material  velocity,  ejected  material  will  simply  fall  back 
to  the  asteroid  or  remain  in  orbit  around  it,  resulting  in  the  object  resuming  its  original 
velocity.  To  approximate  the  maximum  radius  and  mass  that  could  be  deflected,  Eqn.  (5) 
solves  for  Racnt  with  v®  =  0.  At  this  radius,  all  attempts  to  deflect  the  object  are 
impossible  by  this  method  of  explosion.  The  critical  radius  is  defined  as 


R. 


(6) 


where  is  the  average  density  of  the  object  and 


m„ 


(7) 


Finally,  the  critical  mass  is  given  by 


-3 

Vm 
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Knowles  and  Erode  have  published  data  showing  that  a  conservative  estimate  of  the 
ejection  velocity,  using  currently  available  technology,  is  approximately  100  m/s.  Using 
and  average  density  of  5  g  cm'^  for  the  assailant  object,  an  average  ejection  velocity  of 
100  m  s'^  for  current  technology,  the  critical  mass  and  radius  for  the  object  are 
approximately  5  x  10**  kg  and  60  km,  respectively.  Asteroids  with  masses  or  radii 
greater  than  these  values  cannot  be  deflected  with  current  technology,  regardless  of  the 
warning  time  involved.  [Ref.  9] 

2.  Directional  Considerations 

The  orbit  perturbation  requirements  to  deflect  objects  from  collision  with  the 
Earth  are  an  important  consideration  when  planning  a  deflection  mission.  Two  scenarios 
are  considered:  (1)  A  short  time-scale  deflection  where  equations  of  motion  are 
linearized,  and  (2)  A  long  time-scale  deflection,  using  orbital  equations  of  motion.  These 
two  scenarios  are  developed  fully  below,  the  results  of  which  will  be  used  to  generate  the 
change  in  velocity  requirements  for  a  hypothetical  asteroid. 

a.  Short  Time-Scale  Deflection 

The  short  time-scale  deflection  scenario  is  easily  developed  using 
rectilinear  equations  of  motion.  This  is  valid  if  the  time  scale  involved  is  short  compared 
to  the  orbital  period,  P.  To  prevent  a  collision,  a  change  in  velocity  must  be  imparted 
orthogonal  to  the  flight  path  of  the  object.  The  displacement,  6,  that  must  be  achieved  by 
a  change  in  velocity,  Av,  is 

5  =  Av  •  t.  (9) 

Assuming  the  object  were  on  a  path  such  as  to  intersect  the  center  of  mass 
of  the  Earth,  the  object  must  be  deflected  a  minimum  of  one  Earth  radius.  This  provides 
no  margin  of  safety,  however,  and  therefore,  taking  the  minimum  acceptable  distance  to 
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be  two  Earth  radii,  the  minimum  change  in  velocity  required  to  deflect  a  potential  threat 
is  given  by 

2IL  147.6^5-'  (10) 

Av  W - W - - . 

t  t,  days 

b.  Long  Time-Scale  Deflection 

For  a  deflection  to  be  accomplished  over  a  longer  period  of  time,  the 
equations  of  motion  cannot  be  simplified  to  the  rectilinear  case  and  the  two-body  orbital 
mechanics  equations  of  motion  apply.  There  are  three  directions  with  respect  to  the  flight 
path  which  will  cause  an  orbiting  body  to  deviate  from  the  original  orbit:  (1)  orthogonal 
to  the  flight  path,  along  the  line  parallel  to  the  angular  momentum  vector,  (2)  orthogonal 
to  the  flight  path,  in  the  plane  of  the  orbit,  and  (3)  along  the  flight  path,  in  the  orbit  plane. 
A  two  body  circular  orbit  is  assumed  to  simplify  the  problem. 

Given  a  longer  response  time,  the  best  deflection  direction  may  not  be 
orthogonal  to  the  flight  path  direction,  as  in  the  short  time-scale  deflection.  Each  option 
is  assessed  below  to  determine  the  optimal  deflection  direction. 

(1)  Out  of  plane,  orthogonal  to  flight  path.  A  change  in  velocity 
made  orthogonal  to  the  orbit  plane,  as  shown  in  Figure  10, 
results  in  a  change  of  inclination  according  to 


Figure  10.  Deflection  Out  of  Plane,  Orthogonal  to  Flight  Path. 


where  v  is  the  orbital  speed  of  the  object  and  i  is  the  required  inclination  change.  The 
maximum  displacements  occur  n/2  and  Sn/l  around  the  orbit.  The  orbital  period  remains 
imchanged.  Using  a  small  angle  approximation  gives 


Av  =  vi  =  v  — 
r 


where  5  is  the  amount  the  body’s  orbit  will  be  deflected  idl  past  the  point  of  the 
deflection,  and  r  is  the  radius  of  the  object’s  orbit  around  the  Sun.  Finally,  substituting  in 
for  the  orbital  velocity  in  terms  of  orbital  period,  P,  gives 


The  minimum  change  in  velocity  required  to  deflect  the  asteroid  two  Earth  radii  is  given 
by 
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(15) 


^v  = 


P 


For  a  circular  orbit  at  1  AU,  the  minimum  change  in  velocity  required  to  achieve  the 
required  deflection  is  approximately  2.54  m/s  and  should  be  imparted  90°  prior  to  the 
expected  impact  point  in  the  asteroid’s  orbit. 

(2)  In  plane,  orthogonal  to  flight  path.  The  change  in  velocity 
required  to  deflect  an  object  using  a  deflection  in  the  asteroid’s  orbital  plane,  orthogonal 
to  the  flight  path  is  found  in  a  manner  similar  to  (1),  above.  As  before,  the  orbital  period 
remains  the  same.  The  asteroid  is  displaced  along  it’s  track  with  the  maximum 
displacement  occurring  nil  and  37t/2  around  the  orbit,  as  shown  in  Figure  11.  The 
maximum  deflection 


Figure  11.  Deflection  In  Plane,  Orthogonal  to  Flight  Path. 


is 

^2AvP  (16) 

^max  —  _ 

7t 

and  the  minimum  change  in  velocity  required  to  deflect  the  asteroid  two  Earth  radii  is 
given  by 
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(17) 


Avs 


Rji 

P  ■ 


The  required  change  in  velocity  for  this  case  is  approximately  0.635  m  s''. 

(3)  In  plane,  along  flight  path.  A  deflection  along  the  flight  path, 

shown  in 


will  result  in  a  net  increase  or  decrease  in  the  orbital  speed.  If  there  is  a  net  increase  in 
speed,  the  point  at  which  the  deflection  occurred  will  be  the  perihelion  of  a  new  orbit. 
Conversely,  if  there  is  a  net  decrease  in  speed,  the  deflection  point  will  be  the  aphelion  of 
a  new  orbit.  Irrespective  of  deflection  direction,  the  period  of  the  deflected  object’s  orbit 
will  change.  This  change  in  period  serves  as  the  mechanism  by  which,  over  time,  an 
object  can  be  influenced  so  as  to  miss  a  collision  with  the  Earth.  The  amount  the  asteroid 
is  displaced  during  each  orbit  is  directly  proportional  to  the  change  in  period  of  effected 
by  the  deflection.  This  is  given  by 

^per.orbi,=VlSt.  (18) 

The  approximate  displacement  for  a  deflection  along  the  flight  path  can  be  derived  from 
the  equations  for  the  specific  mechanical  energy  of  a  body, 

P  P  -fx  (19) 

the  orbital  period  of  a  body. 
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and  the  circular  velocity  for  an  orbit, 


Taking  the  derivative  of  Eqn  (19)  and  substituting  for  p/a  gives 

Av  Aa 


while  taking  the  derivative  of  Eqn.  (20)  and  resubstituting  in  for  T  gives 


AT  =  3T^. 
2a 


(21) 


(22) 


(23) 


Substituting  (22)  and  (23)  into  (18)  gives  the  maximum  deflection  per  orbit,  in  terms  of 
the  original  orbital  period  of  the  asteroid  and  the  change  in  velocity  imparted  to  the 
object,  as 

V-./,=3rAv.  (24) 


So,  the  change  in  velocity  required  to  deflect  an  object  two  Earth  radii  is  approximately 


where  n  is  the  number  of  orbits  prior  to  impact.  If  n  =  1,  the  change  in  velocity  is 
approximately  0.135  m  s'*.  This  is  better  than  either  case  (1)  or  (2).  If  the  impact  can  be 
predicted  as  much  as  a  decade  in  advance,  the  change  in  velocity  required  to  deflect  the 
asteroid  is  reduced  to  1.35  cm  s‘^  This  is  an  order  of  magnitude  less  than  case  (2)  and 
two  orders  of  magnitude  less  than  case  (1)  and  is  an  achievable  quantity  in  terms  of 
technological  feasibility. 
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IV.  DEVELOPMENT  OF  HYPOTHETICAL  ORBITS 


A.  INTRODUCTION 

In  order  to  accurately  test  the  sensitivity  of  the  change  in  velocity  requirements 
presented  in  section  III,  accurate  orbital  parameters  for  hypothetical  asteroid  objects  were 
first  developed.  This  section  first  describes  the  coordinate  systems  used  and  the 
characteristics  of  circular,  elliptical,  parabolic,  and  hyperbolic  orbits.  Next,  this  section 
describes  the  process  by  which  the  Earth’s  position  in  heliocentric  coordinates,  at  any 
given  time,  are  used  to  generate  the  heliocentric  coordinates  of  an  intersecting  asteroid 
orbit.  The  size,  shape,  and  inclination  of  the  asteroid  orbit  are  variables  defined  by  the 
user  with  the  only  limitation  being  that  the  two  orbits  intersect.  MATLAB  scripts  which 
generate  these  procedures  are  contained  in  the  Appendix.  Orbits  can  be  generated  for 
circular,  elliptical,  parabolic,  or  hyperbolic  trajectories. 

Development  of  hypothetical  orbits  followed  a  rigorous  routine  allowing  for 
repeatability  of  the  work  done  while  minimizing  the  workload  for  completing  multiple 
simulations.  The  methodology  and  theory  used  in  the  development  of  the  MATLAB 
scripts  are  discussed  below. 

B.  COORDINATE  SYSTEMS 

The  first  requirement  for  describing  an  orbit  is  to  define  a  suitable  reference 
fi'ame.  In  many  cases  this  means  finding  an  appropriate  inertial  coordinate  system.  For 
orbits  around  the  Sun  such  as  planets,  asteroids,  and  comets,  the  heliocentric-ecliptic 
coordinate  system  is  often  used.  Satellites  in  orbit  around  the  Earth  use  the  geocentric- 
equatorial  system.  [Ref  1 0]  Another  convenient  coordinate  frame  used  for  describing  the 
orbit  of  a  body  is  the  perifocal  coordinate  system.  Each  rectangular  coordinate  frame  is 
defined  by  specifying  the  origin,  the  fundamental  plane  (i.e.  the  X-Y  plane),  and  the 
direction  of  each  axis.  In  developing  the  hypothetical  orbits  used  for  the  simulations,  the 
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heliocentric-ecliptic  and  perifocal  coordinate  systems  were  used  extensively.  These  are 
described  in  detail  below. 

1.  Heliocentric  Ecliptic  Coordinate  System 

The  heliocentric-ecliptic  coordinate  system,  as  the  name  implies,  has  its  origin  at 
the  center  of  the  Sun.  The  letters  XYZ  describe  the  three  principle  axes  as  shown  in 
Figure  13.  The  fundamental  plane,  the  X-Y  plane,  coincides  with  the  ecliptic,  which  is 
the 


Figure  13.  The  Heliocentric-Ecliptic  Coordinate  System.  From  Ref  10. 

plane  defined  by  the  Earth’s  orbit  around  the  Sun.  The  line  of  intersection  of  the  Earth’s 
equatorial  plane  and  the  ecliptic  plane  (XY)  where  the  Sim  crosses  the  equator  from  south 
to  north  in  its  apparent  annual  motion  along  the  ecliptic  is  known  as  the  vernal  equinox. 
When  the  vector  from  the  Earth  to  the  Sun  (X-axis)  points  towards  what  is  now 
approximately  the  constellation  of  Pisces,  it  coincides  with  the  vernal  equinox,  also 
known  as  the  first  day  of  spring.  The  Y  axis  is  90°  from  the  Z  axis  in  the  direction  of  the 
Earth’s  rotation  around  the  Sun.  The  Z  axis  completes  the  right-handed  coordinate 
system.  [Ref.  10] 

2.  Perifocal  Coordinate  System 

The  perifocal  coordinate  system,  shown  in.  Figure  14,  is  one  of  the  most  useful 
coordinate  systems  for  describing  the  motion  of  a  body.  The  letters  PQW  are  used  to 
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describe  the  three  principle  axes.  The  fundamental  plane,  the  PQ  plane,  is  in  the  plane  of 
the  body’s  orbit  and  the  origin  is  at  the  focus  of  the  primary  body,  in  our  case  the  Sun. 
The  P  axis  is  in  the  direction  of  perihelion  while  the  Q  axis  is  90°  from  the  P  axis  in  the 
direction  of  rotation.  The  W  axis  completes  the  right-handed  coordinate  system. 


Figure  14.  Perifocal  Coordinate  System. 


C.  CLASSICAL  ORBITAL  ELEMENTS 

Five  classical  orbital  elements,  also  known  as  Keplerian  elements,  are  sufficient  to 
describe  the  size,  shape,  and  orientation  of  an  orbit  A  sixth  element  is  required  to  locate 
a  body  in  that  orbit. 

There  are  many  ways  of  locating  the  position  of  a  body  in  the  orbit  plane  which 
may  be  substituted  for  time  of  perihelion  passage. 

1.  True  anomaly  at  epoch,  Vq,  is  the  angle  in  the  plane  of  the  body’s  motion 
measured  from  perihelion  to  the  position  to  the  body  at  a  given  time  (epoch). 

2.  Argument  of  latitude  at  epoch,  Uo,  is  the  angle  in  the  plane  of  the  orbit  between 
the  longitude  of  the  ascending  node  and  the  radius  vector  to  the  body  at  a  given  time.  If  co 
and  Vo  are  both  defined  then 

Uo=(0+Vo.  (26) 
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If  there  is  no  ascending  node,  as  in  an  equatorial  orbit,  then  both  (o  and  Uo  are 
undefined. 

3.  True  longitude  at  epoch,  lo,  is  the  angle  measured  between  X  and  r,  the  radius 
vector  to  the  body,  measured  eastward  to  the  ascending  node,  if  it  exists,  and  then  in  the 
orbital  plane  to  r.  If  Q,  co,  and  Vo  are  all  defined,  then 

lo=Q+G)+Vo=n+Vo=Q+Uo.  (27) 

Two  other  terms  used  in  the  descriptions  of  orbits  are  direct,  or  prograde,  and 
retrograde.  Direct  means  easterly,  in  accordance  with  the  right-hand  rule,  and  is  the  same 
direction  in  which  the  Sun,  Earth,  and  most  of  the  planets  and  their  moons  rotate  on  their 
axes  and  the  direction  in  which  all  planets  rotate  around  the  Sun.  A  retrograde  orbit,  as 
the  term  implies,  is  the  opposite  of  direct.  Table  4  provides  a  summary  of  the  various 
types  of  orbits  and  the  orbital  parameters  that  result. 


Parameter 

Circle 

Parabola 

Eccentricity 

e=0 

0<e<l 

e=l 

e>l 

Energy 

-fi/2a<0 

E=0 

mu/2a>0 

Period 

P=27i:/n 

P=27t/n 

- 

- 

Anomaly 

Undefined  or 
arbitrary 

Eccentric: 

^  V  fl  +  eV  E 

tan-=  - —  tan— 

2  vl  —  e/  2 

Parabolic,  D 
v  D 

2  ^ 

Hyperbolic,  F 

.  V  fl  +  eV 

Table  4.  Keplerian  Orbits 


D.  DERIVATION  OF  FORMULAE 
1.  Procedure 

a.  Orbital  Elements  From  r  and  v 

Conversion  from  position  and  velocity  vectors  to  orbital  elements  is  done 
using  standard  procedures.  For  the  purposes  of  this  example,  the  Earth’s  orbital 
parameters  are  being  determined  from  position  and  velocity  vectors. 
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First  the  angular  momentum  is  determined.  From  the  relative  form  of  the 


basic  two-body  equation 


Cross  multiplying  both  sides  by  the  position  vector,  r,  gives 


-  -  -  F  - 
rxr  +rx^r 
r 


Expanding  this  equation,  one  then  finds  the  angular  momentum 

h  =  fxv  =  constant.  (30 

The  vector  pointing  to  the  node  is  then  calculated  using 

n  =  Kxh  (31 

where  K  is  the  unit  vector  in  the  Z  direction.  If  the  magnitude  of  the  node  vector  is  zero, 
the  orbit  is  in  the  fundamental  plane.  The  MATLAB  script  used  to  calculate  this  vector 
has  a  filter  to  set  the  vector  equal  to  zero  if  the  magnitude  is  less  than  1x10'^. 

Next  the  eccentricity  vector  is  calculated  with 


e=—  v^-—  F-(F-v)v  . 

p  1_V  rJ 


The  magnitude  of  this  vector  provides  the  eccentricity.  It  is  important  to  realize  that  this 
vector  is  zero  for  circular  orbits.  o 

The  semi-latus  rectum,/?,  is  calculated 


as  are  the  semi-major  axis 


and  the  radius  at  aphelion  and  perihelion 


K 


(35) 


-A., 

\-e’  ”  l  +  e' 


Once  the  size  and  shape  parameters  have  been  determined,  all  that  remains 
is  to  find  the  orientation  of  the  Earth’s  orbital  plane  and  the  position  of  the  Earth  in  that 
orbit.  First,  the  inclination  is  found  with 


cos(/)  = 


h-K 

\h\ 


(36) 


There  is  no  need  to  check  for  angle  ambiguities  as  the  inverse  cosine  always  returns  an 
angle  between  0°  and  1 80°. 


Next,  the  longitude  of  the  ascending  node,  which  varies  from  0°  to  360° 
must  be  determined.  By  definition,  the  Earth  is  in  the  mean  ecliptic  and  the  inclination  is 
zero  deg.  By  convention,  the  longitude  of  ascending  node  is  defined  to  be  zero  deg;  there 
is  no  real  “ascending  node”.  However,  due  to  small  perturbations  in  the  Earth’s  orbit 
around  the  Sun,  the  center  of  mass  of  the  Earth  deviates  out  of  the  mean  ecliptic  to  what 
is  known  as  the  true  ecliptic.  As  a  result,  when  there  are  instantaneous  components  of  the 
position  and  velocity  vector  normal  to  the  mean  ecliptic,  very  small,  but  calculable  values 
for  inclination  can  be  found.  Also,  a  value  for  the  longitude  of  ascending  node  can  be 
found.  If  the  instantaneous  position  and  velocity  vectors  lie  in  the  plane  of  the  ecliptic, 
the  longitude  of  ascending  node  and  inclination  are  zero  deg.  and  the  only  rotation  is  that 
of  the  argument  of  perihelion.  As  mentioned  previously,  orbits  aligned  with  the 
fimdamental  plane,  in  our  case  the  ecliptic,  have  no  node  vector  and  the  longitude  of  the 
ascending  node  is  undefined.  However,  a  value  for  the  majority  of  orbits  is  obtained  from 


cos(Q)  = 


I  ■  n 

|fNi 


(37) 


In  this  case  a  quadrant  check  must  be  conducted  and  if  the  dot  product  of  J  and  n  is  less 
than  zero,  then  Q  must  lie  between  180°  and  360°. 
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The  argument  of  perihelion,  which  also  varies  from  0°  to  360°  is 
calculated.  Again,  this  value  is  undefined  for  circular  orbits  and  orbits  in  the  ecliptic 
because  perihelion  does  not  exist  for  the  former  and  there  is  no  node  vector  in  the  later. 


(38) 


Angle  ambiguities  are  resolved  by  computing  the  dot  product  of  K  and  e.  If  this  value  is 
less  than  zero,  then  perihelion  is  below  the  ecliptic  and  the  argument  of  perihelion  is 
between  1 80°  and  360°. 

The  final  classical  orbital  element  is  the  true  anomaly  at  epoch,  which  can 
vary  from  0°  to  360°.  As  mentioned  previously,  the  true  anomaly  is  the  angle  in  the  orbit 
plane  measured  from  perihelion  to  the  Earth’s  current  position.  The  true  anomaly  is  then 


(39) 


The  correct  quadrant  for  the  true  anomaly  is  foimd  from  the  dot  products  of  r  and  v.  If 
the  dot  product  is  less  than  zero,  then  the  true  anomaly  is  between  180°  and  360°. 


b.  Singularity  Solutions 

Several  special  cases  must  be  considered,  especially  when  calculating 
values  for  the  orbital  elements  using  a  computer  routine.  The  first  occurs  when  die 
orbital  plane  is  the  same  as  the  fundamental  plane.  This  is  especially  relevant  to  the  orbit 
of  the  Earth  aroimd  the  Sun.  By  definition,  the  Earth  is  in  the  ecliptic,  the  inclination  and 
longitude  of  ascending  node  are  zero  and  the  normal  equation  to  calculate  the  argument 
of  perihelion  cannot  be  used.  However,  the  instantaneous  position  and  velocity  vectors 
may  contain  small  components  out  of  the  plane  of  the  mean  ecliptic.  There  are  small 
values  for  inclination  and  longitude  of  the  ascending  node.  To  obtain  the  required 
accuracy,  these  angles  are  included  when  calculating  the  asteroid’s  orbit.  Appropriate 
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filters  handle  the  case  when  the  Earth’s  position  lies  exactly  in  the  mean  ecliptic,  and  the 
inclination  and  longitude  of  ascending  node  are  zero. 

These  singularities  are  removed  by  using  the  true  longitude  of  perihelion, 
TiJtnie,  which  Combines  Q.  and  co.  It  is  found  through 

,  he  (40) 

w 

which  is  valid  for  elliptical  orbits  which  lie  in  the  ecliptic.  Again,  a  quadrant  check  is 
necessary  and  is  found  using  the  dot  product  of  J  and  e.  If  this  value  is  less  than  zero, 
the  true  longitude  of  perihelion  is  between  180°  and  360°. 

c.  Orbital  Anomaly  Determination 

There  are  several  other  useful  formulas  that  must  be  derived  prior  to 
writing  script  files  for  the  plotting  of  orbits  and  calculation  of  orbital  parameters.  First, 
for  elliptical  orbits,  such  as  those  of  the  Earth,  EGAs,  and  short-period  comets,  formulas 
for  determining  the  eccentric  anomaly  and  the  true  anomaly  prove  quite  useful.  Figure  15 
shows  an  x,y  Cartesian  coordinate  system  centered  at  C.  F  is  the  focus  of  an  ellipse,  in 
this  case  the  Srm.  An  auxiliary  circle  of  radius  a  with  center  C  is  constructed.  P  is  the 
location  of  the  Earth  in  its  orbit  about  the  Sun  and  Q  is  the  point  where  the  perpendicular 
to  the  semi-major  axis  intersects  the  auxiliary  circle.  The  angle  E  is  the  eccentric 
anomaly  and  v  is  the  true  anomaly  of  the  point  P. 


44 


Figure  1 5.  Orbital  Anomalies  for  Elliptic  Motion 


Another  identity.  Gauss’  law,  provides  a  useful  relation  between  v  and  E. 
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1 

tan— V  = 
2 


(47) 


l  +  e  1 

- - tan-E. 

1-e  2 


Because  1/2  v  and  1/2  E  are  always  in  the  same  quadrant,  there  is  no  angle  ambiguity  to 
resolve  when  solving  this  form  of  the  equation. 


2.  Determination  of  r  and  v  From  Intersecting  Orbit  Parameters 

From  a  given  orbit  and  the  position  of  a  body  within  that  orbit,  a  second  orbit  can 
be  constructed  which  intersects  the  first  orbit.  The  location  of  a  body  within  the  second 
orbit  which  intersects  the  body  in  the  first  orbit  is  specified  by  the  formulation  of  the 
solution.  Values  for  semi-major  axis,  eccentricity,  and  inclination  are  variables  and  can 
be  any  value  with  the  only  limitation  that  the  two  orbits  must  have  at  least  one  point  of 
intersection.  Finally,  the  position  and  velocity  vectors  in  the  new  orbit  at  the  point  of 
intersection  are  calculated.  These  position  and  velocity  vectors  are  in  perifocal 
coordinates.  The  detailed  procedure  for  calculating  these  vectors  is  presented  below. 

a.  Circular,  Elliptical,  and  Hyperbolic  Orbits. 

The  polar  equation  for  an  ellipse 

a(l-e')  (48) 

r  = - 

1  +  ecosv 

can  be  used  to  solve  for  the  true  anomaly  at  epoch.  The  magnitude  of  r  in  this  equation  is 
the  magnitude  of  the  Earth’s  position  vector  at  epoch.  Once  the  true  anomaly  at  epoch  is 
known,  the  position  vector  in  the  perifocal  coordinate  system  can  be  determined  using 

F  =  rcosvP  +  rsinvQ  (49) 

where  the  magnitude  of  r  is  once  again  the  magnitude  of  the  Earth’s  position  vector  and 
epoch. 

Differentiating  Eqn.  (49)  in  this  “inertial”  perifocal  frame  gives 
V  =  r  =  (j  cosv  -  rv  sinv)P  +  (r  sinv  -t-  rv  cosv)Q .  (50) 


46 


Substituting  in 


gives 


/  * 

rv  =  J— (1  +  ecosv) 
P 


V  =  I— [-sinvP  +  (e+  cosv)Q] 
P 


(51) 


(52) 


b.  Parabolic  Trajectories. 

The  semi-major  axis  of  a  parabolic  orbit  is  undefined.  If  the  desired  orbit 
is  parabolic  in  shape,  the  parameter,/?,  is  substituted  for  the  semi-major  axis.  The  polar 
equation  of  the  parabola  is  then 

r  =  -P—. 

1  -i-cosv  ’ 

This  equation  is  solved  for  v,  as  discussed  previously.  The  formulas  for 
calculating  the  position  and  velocity  vectors  in  perifocal  coordinates  are  then  identical  to 
the  elliptical  case  discussed  previously. 


3.  Orbit  Rotations 

Rotations  between  coordinate  frames  are  accomplished  using  Direction  Cosine 
Matrices  (DCMs).  Coordinate  rotations  may  be  accomplished  by  any  number  of 
methods.  One  common  method  of  conducting  coordinate  rotations  is  through  a  3-1-3 
rotation:  A  rotation  about  the  3-axis,  followed  by  a  rotation  about  the  new  1-axis,  and 
finally  a  3-rotation  about  the  new  3-axis.  And  example,  shown  in  Figure  16,  converts 
fi-om  perifocal  coordinates  to  Heliocentric-ecliptic  coordinates.  The  rotation  matrices  as  a 
function  of  angle  of  rotation,  a,  are  presented  below.  The  subscript  associated  with  each 
transformation  is  the  axis  about  which  the  rotation  is  to  be  accomplished. 
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(54) 


Figure  16.  Coordinate  Transformations.  From  Ref.  10. 
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V.  CASE  STUDIES. 


A.  INTRODUCTION 

Two  tests  were  conducted  to  analyze  the  effects  of  orbital  perturbations  on  an 
assailant  asteroid.  The  magnitude  of  the  velocity  applied,  is  the  change  required  along  the 
flight  path,  to  change  the  orbital  period  of  a  circular  orbit  such  that  the  change  in  position 
after  one  orbit  is  equal  to  two  Earth  radii.  This  magnitude  is  approximately  13.5 
cm/sec/orbit.  The  first  used  a  solar  system  simulation  modeling  program  called  Dance  of 
the  Planets.  [Ref  4]  This  program,  described  in  further  detail  below,  allows  the  user  to 
input  hypothetical  objects  into  a  relatively  high  fidelity  simulation  of  the  solar  system.  It 
includes  gravitational  models  for  the  Sim  and  all  of  the  planets.  Objects  can  be  imported 
into  the  program  using  their  orbital  elements  or  their  position  and  velocity  vectors  at  a 
given  time.  The  model  then  calculates  new  orbital  positions  and  velocities  for  each  body 
included  in  the  model  using  a  full  gravitational  model.  The  simulation  accurately 
computes  both  forward  and  backward  in  time.  Although  the  program  provides  a 
reasonable  model  visualization,  an  inaccuracy  was  discovered  during  testing  which 
prevented  importing  of  position  and  velocity  vectors  of  sufficient  accuracy  for  valid  tests 
to  be  conducted.  This  inaccuracy  is  demonstrated  and  explained  in  section  7  of  this 
chapter.  The  discussion  focuses  on  the  procedure  used  to  generate  the  desired  asteroid 
vectors.  This  would  prove  useful  for  follow-on  research  should  the  program  be  modified 
to  provide  the  required  accuracy. 

The  second  method  used  to  evaluate  sensitivity  of  deflection  distance  and 
direction  to  variations  in  orbit  eccentricity  and  size  (semi-major  axis)  consisted  of  a 
relatively  simple  model  to  determine  miss  distance  based  on  the  asteroid’s  mean  anomaly, 
M,  and  mean  motion,  n.  The  position  and  velocity  vectors  at  some  time  prior  to  collision 
were  calculated,  new  orbital  elements  determined,  and  the  new  position  at  the  original 


49 


time  found.  The  distance  between  the  center  of  the  Earth’s  orbit  and  the  asteroid’s  orbit 
was  then  calculated.  The  effect  of  varying  specific  orbital  parameters  while  keeping  the 
magnitude  of  the  change  in  velocity  constant  were  analyzed.  Parameters  varied  were 
orbital  eccentricity,  semi-major  axis,  and  deflection  direction. 

B.  DANCE  OF  THE  PLANETS 

1.  Hypothetical  Orbit  Determination 

The  purpose  of  this  section  is  to  present  the  process  used  to  calculate  asteroid 
object  orbital  elements  and  heliocentric  coordinates  from  the  Earth’s  heliocentric 
coordinates.  The  position  and  velocity  vectors  obtained  are  the  position  of  the  asteroid  in 
its  new  orbit  around  the  sun  in  heliocentric  coordinates.  This  is  exactly  the  same  position 
as  that  of  the  Earth,  hence  collision.  The  asteroids  velocity  vector  is  that  vector  which 
meets  the  size  and  shape  requirements  dictated  earlier  in  the  orbit  development  process 
and  is  the  asteroid’s  instantaneous  velocity  vector  at  collision.  These  vectors  are  then 
imported  into  Dance  of  the  Planets  at  the  epoch  for  which  the  Eeirth’s  original  r  and  v 
were  obtained.  The  orbit  was  then  propagated  backwards  in  time.  New  velocity  vectors 
with  a  delta-v  applied  were  calculated  and  imported  into  the  program.  These  new  vectors 
were  propagated  forward  to  determine  if  collision  had  been  averted.  The  procedure, 
summarized  here,  is  described  in  detail  in  the  sections  which  follow. 

a.  First,  the  Earth’s  orbital  elements  are  calculated  from  r  and  v. 
Specifically,  true  anomaly,  v,  is  required  so  that  the  position  in  the 
orbit  is  known. 

b.  The  asteroid’s  orbit  size,  shape,  and  inclination  (a,  e,  and  i)  are 
selected. 

c.  The  asteroid’s  orbital  elements  are  calculated  in  perifocal  coordinates. 

d.  The  asteroid’s  orbit  in  perifocal  coordinates  is  rotated,  using  a  3-1-3 
rotation,  to  make  it  coplanar  with  the  Earth’s  orbit.  The  eccentricity 
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vectors  are  aligned.  This  makes  the  Sun  and  the  two  orbit’s 
perihelions  collinear. 

e.  A  second  3-1-3  rotation  is  performed  to  obtain  the  desired  inclination 
and  align  the  two  r  vectors.  The  first  3 -rotation  is  through  an  angle 
Vearth  and  aligns  the  asteroid’s  perihelion  with  the  position  of  the  Earth 
in  its  orbit.  The  1 -rotation  inclines  the  orbit  to  the  desired  value.  The 
final  3-rotation,  through  an  angle  -Vasteroid,  aligns  the  asteroid  with  the 
Earth’s  current  position. 


2.  Initialization. 

The  first  step  in  the  development  of  hypothetical  orbits  and  completing  the 
simulation  of  Near  Earth  Objects  (NEOs)  was  to  select  a  desired  date  and  time  for  the 
impact.  For  the  purposes  of  this  research,  the  date  chosen  was  January  01, 2007  at  00:00 
UT.  This  date  was  selected  as  it  gave  a  real  world  flavor  to  the  problem,  approximately 
one  decade  would  elapse  from  the  time  of  discovery  to  the  time  of  impact.  Universal 
Time  (UT)  is  the  standard  time  in  the  time  zone  centered  on  zero  degrees  longitude  in 
Greenwich,  England.  UT  is  used  for  calculating  where  an  object  is  in  the  sky  relative  to 
the  horizon  at  a  particular  location. 

Next  the  position  and  velocity  vectors  for  Earth  at  the  desired  impact  date  were 
obtained  firom  Dance  of  the  Planets.  The  operator’s  manual  contains  an  excellent  tutorial 
and  provides  the  information  required  to  operate  the  program.  Several  preliminary  steps 
must  be  completed  prior  to  date  initialization.  The  simulation  is  initialized  firom  Space 
View.  Second,  from  the  main  access  screen,  any  planets  that  are  to  be  excluded  fi'om  the 
gravitational  model  are  deselected.  This  is  indicated  by  a  zero  in  the  on  column  next  to 
the  planet  name.  Returning  to  the  simulation  Space  View,  confirm  the  simulation  pace  is 
set  to  true  time.  This  minimizes  the  elapsed  time  upon  returning  to  the  simulation  after 
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initializing  the  date.  In  addition,  this  also  minimizes  the  change  in  planetary  position  and 
velocity  vectors  during  the  time  from  initialization  to  reselecting  the  main  access  facilities 
display  where  the  Earth’s  vectors  are  obtained.  Finally,  enter  the  desired  date.  Digital 
hours  and  minutes  appear  at  low  simulation  paces.  Local  time,  determined  by  site 
longitude,  is  shown  unless  time  is  offset. 

Immediately  after  entering  the  desired  date,  reenter  the  main  access  screen.  For 
each  planet  selected,  Dance  of  the  Planets  calculates  the  period,  the  distance  to  the  Earth 
and  Sun  in  astronomical  units  (AU),  heliocentric  longitude  (calculated  and  simulated), 
right  ascension  of  the  ascending  node,  declination,  and  apparent  magnitude.  The 
difference  between  calculated  and  simulated  heliocentric  longitude  is  that  the  calculated 
value  is  the  position  as  determined  from  ephemeris  equations  while  the  simulated  value 
gives  the  position  in  “simulator  space.”  Comparison  of  the  two  values  provides  an 
accuracy  check  for  the  user  between  the  ephemeris  equations  and  the  Dance  of  the 
Planets  model. 

The  instantaneous  position  and  velocity  vectors  of  the  planets  provide  the 
essential  information  to  describe  dynamic  state.  With  these  datum  the  six  orbital 
parameters  of  eccentricity,  semi-major  axis,  right-ascension  of  the  ascending  node, 
inclination,  argument  of  perihelion,  and  true  anomaly,  can  be  calculated. 

3.  Constants 

The  program  provides  the  user  with  the  option  of  obtaining  the  position  and 
velocity  vectors  or  the  orbital  elements  at  a  given  epoch.  The  display  shows  the 
calculated  vectors  for  the  planets  in  heliocentric  coordinates,  XYZ,  as  well  as  the  epoch 
Julian  date.  Units  for  the  position  and  velocity  vectors  are  E6  km  and  km/sec, 
respectively.  The  epoch  date  is  presented  in  the  form  24ddddd.ffff  and  is  reference  to  the 
epoch  (equinox)  of  the  year  2000  (J2000),  as  are  contemporary  star  charts.  Only  Julian 
dates  beginning  with  24ddddd  can  be  used,  giving  a  14,680  year  vwndow  for  simulations. 
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[Ref.  4]  Mean  planetary  constants  for  epoch  J2000  of  the  Earth  are  given  in  Table  5. 
Solar  constants  are  shown  in  Table  6. 


Parameter 

Value 

Semi-major  axis 

AU 

km 

1.000  001  017  8 

149,598,023 

Eccentricity 

0.016  708  617 

Inclination  ~  deg. 

0.000  000  000 

Long,  of  ascending  node  (Q)  ~  deg 

0.000  00000 

Long,  of  perihelion  (ro)  ~  deg 

102.937  348  08 

True  longitude  (L)  ~  deg 

100.466448  51 

Orbital  period  (P)  ~  years 

0.999  978  62 

Orbital  velocity  (v)  ~  km/s 

29.77859 

Equatorial  radius  (Re)  ~  km 

6378.137 

Gravitational  parameter  (p)  ~  km^/s^ 

3.986x10^ 

Mass  (Me)  ~  kg 

5.9742x10^^ 

Rotation  ~  days 

0.997  269  69 

Inc.  of  equator  to  ecliptic  ~  deg 

23.45 

Density  ~  gm/cm^ 

5.515 

Table  5.  Mean  Earth  Constants  for  Epoch  J2000 


Parameter 

Value 

Radius  of  the  Sxm  (Rs)  ~  km 

696,000.000  . 

1.0AU~km 

149,597,870.0 

1 .0  AU/TUsun  ~  km/solar  s 

29.784  691  674  9 

1.327  124  28x10" 

1 .0  TU  ~  solar  days 

58.132  440  9 

Table  6.  Solar  Constants 


4.  Orbital  Element  Selection 

The  orbital  elements  selected  for  simulation,  shown  in  Table  7,  were  selected 
from  a  list  of  all  known  ECA’s  given  by  Rabinowitz  et.  al.  [Ref  5]  and  [Ref  7].  These 
were  selected  because  they  each  have  a  point  of  close  approach  distance  within  20  lunar 
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radii  over  the  course  of  the  next  century,  and  because  they  are  representative  of  each  class 
of  asteroid  as  well  as  a  short-period  comet. 


No.  and  Name 

H 

Approx, 
diamete 
r  (km) 

Depth 

(AU) 

q 

(AU) 

(AU) 

e* 

(deg) 

2340  Hathor 

20.26 

0.2 

0.356 

0.464 

0.844 

5.85 

4660  Nereus 

18.30 

1 

0.092 

0.954 

1.490 

1.43 

4179  Toutatis 

14.0 

5 

- 

0.542 

2.505 

0.640 

0.47 

Encke 

9.8 

- 

!  - 

0.331 

2.283 

0.855 

11.9 

Note:  (1)  Values  used  in  the  development  o: 


f  the  hypothetical  orbit 


Table  7.  Selected  Near-Earth  Objects  Used  for  Simulation  Study. 


5.  Asteroid  Orbit  Rotations 

Two  series  of  3-1-3  rotations  are  used  to  achieve  the  desired  orbit.  The  first  series 
of  rotations  transforms  the  position  and  velocity  vectors  from  perifocal  coordinates  to  the 
Earth’s  orbital  plane.  These  angles  were  previously  calculated  and  are  the  Q,  i,  co  rotation 
angles.  Because  the  true  ecliptic  is  used,  the  Earth’s  orbit  with  respect  to  the  mean 
ecliptic  is  foimd  using  the  three  rotations  calculated  in  Eqns.  (37),  (38),  and  (39), 
inclination,  longitude  of  ascending  node,  and  argument  of  perihelion,  respectively. 

Once  the  assailant  object’s  orbit  is  coplanar  with  that  of  Earth,  and  the  eccentricity 
vectors  are  aligned,  a  second  3-1-3  coordinate  transformation  is  accomplished.  This 
series  of  rotations  serves  two  important  purposes:  (1)  they  rotate  the  orbit  to  the  desired 
inclination,  and  (2)  they  align  the  position  vectors  of  the  Earth  and  the  asteroid. 

The  first  3 -rotation,  through  the  angle  Veanh,  aligns  the  object’s  eccentricity  vector, 
which  points  towards  perihelion,  with  the  Earth’s  position  vector.  The  l-rotation,  about 
the  eccentricity  vector,  effects  the  desired  value  of  inclination.  The  final  3 -rotation  about 
the  orbit  normal,  though  an  angle  -Vasteroid,  aligns  the  asteroid’s  position  vector  with  the 
Earth’s  position  vector. 
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6.  Simulation 


Once  the  position  and  velocity  vectors  are  calculated,  they  are  inserted  into  Dance 
of  the  Planets  as  object  vectors.  The  fictitious  orbit,  at  epoch,  is  one  having  the  same  r  as 
Earth  but  with  a  velocity  vector  as  provided  by  MATLAB  to  provide  the  desired  orbit  size 
and  shape. 

The  simulation  is  run  backwards  in  time  to  an  arbitrary  time.  This  will  be  the  new 
start  time  of  the  simulation.  At  this  point,  the  simulation  is  paused,  and  the  new  orbital 
elements  are  calculated  for  the  asteroids  current  orbital  set.  These  are  not  the  same  as  the 
original  orbit’s  parameters  due  to  the  influence  of  third  bodies  on  the  asteroid’s  orbit 
during  the  simulation  period.  Specifically,  the  influence  of  the  Earth  on  the  asteroid’s 
orbit  during  the  first  month  of  the  reverse  simulation  has  a  significant  impact  on  the  orbital 
parameters. 

Next,  another  MATLAB  script  file  uses  the  current  orbital  parameters  to 
recalculate  the  velocity  vector  of  the  assailant  body  for  a  desired  deita-V  applied  to  the 
body  at  a  given  time.  This  procedure  is  described  in  detail  below. 

Finally,  the  simulation  was  run  forward  in  time  to  determine  if  Dance  of  the 
Planets  achieved  the  simulation  accuracy  required  to  show  the  deflection  of  the  asteroid 
and  determine  if  the  calculated  deflection  is  sufiBcient  to  avoid  collision.  The  results  are 
presented  in  the  next  section. 

All  MATLAB  scripts,  presented  in  the  Appendix,  were  written  in  modular  fashion. 
Running  the  script  from  the  main  program  accesses  additional  script  files  as  required. 

Each  file  can  also  be  run  from  the  Command  window  by  typing  the  function  name  with  the 
appropriate  arguments  provided.  Help  files  are  provided  for  each  function  file.  The 
MATLAB  scripts  also  plot  the  Earth  and  asteroid  orbits.  Aphelion,  perihelion,  the 
position  of  each  body  in  their  orbits  is  indicated  on  the  plots. 
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7.  Results  and  Discussion. 


Following  the  development  of  the  MATLAB  function  files  and  orbital 
development  routines,  a  problem  was  discovered  with  the  accuracy  of  files  imported  to 
Dance  of  the  Planets.  Data  files  used  to  import  elements  into  the  simulation  are  space 
delimited  and  of  the  form  : 

Name  x  y  z  Vx  Vy  Vz  epoch 
The  position  coordinates  were  in  units  of  E6  km  and  contained  up  to  5  decimal  place 
precision.  Velocity  vector  components  were  in  km/sec  and  contained  up  to  6  decimal 
place  precision.  Epoch  date  was  in  Julian  days  and  allowed  inputs  up  to  four  decimal 
places. 

Importing  a  data  file  with  this  level  of  exactness  would  have  been  sufficient  to 
accurately  test  the  deflection  hypothesis  presented  in  section  III.  However,  when  data 
was  imported,  the  values  obtained  were  slightly  different  from  the  original  file.  Table  8 
show  the  result  of  importing  the  data  file  for  asteroid  2340  Hathor.  The  file  contained  the 
consisted  of  the  position  and  velocity  vectors  shown  under  Hathor,  below,  and  was 
obtained  by  exporting  actual  position  and  velocity  data  from  Dance  of  the  Planets.  The 
example  column  listed  next  to  Hathor  shows  how  the  values  differed  following 
importation. 
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Ast. 


Hathor 


Example 


X  U292Q1  13.29217 

y  -75.41408  -75.41404 

z  0.02352  0.02351 

Vx  44.735765  44.735756 

Vy  19.744673  19.744721 

Vz  -4.868622  -4.868622 

Table  8.  Imported  Asteroid’s  Altered  Position  and  Velocity  Vectors. 

The  magnitude  of  the  change  between  the  desired  values  and  the  values  as  altered 
by  Dance  of  the  Planets  is  shown  in  Table  9.  If  the  values  had  been  imported  correctly 
the  position  and  velocity  vectors  would  be  identical. 


Asteroid 

Delta  r  (km) 

Delta  V  (m/s) 

108.1665 

0.1082 

Table  9.  Position  and  Velocity  Change  from  Original  Orbit. 

While  these  errors  were  relatively  small,  they  were  of  sufficient  magnitude  to 
impart  significant  errors  at  the  initiation  of  the  simulation  and  prevented  precise  and 
accurate  testing  of  the  theories;  the  magnitude  of  the  change  in  velocity  to  be  imparted  to 
the  asteroid  was  of  approximately  the  same  magnitude  as  the  error  generated  by  the 
program.  Additionally,  these  errors  were  not  consistent  in  magnitude  or  direction  from 
one  test  to  the  next.  The  errors  could  not  be  biased  out  by  altering  the  date  of  epoch  for 
the  simulation. 

Discussions  with  Terry  Hannon  of  Applied  Research  &  Consulting,  Inc.,  one  of 
the  software  developers  for  the  program,  lead  to  speculation  that  the  problem  may  lie  in 
the  precision  of  the  epoch  currency.  The  program  presents  accuracies  to  one  ten- 


57 


thousandth  of  a  day,  approximately  eight  seconds,  while  the  accuracies  required  for  the 
thesis  testing  was  on  the  order  of  one  second  accuracy.  In  fact,  although  the  epoch  date  is 
displayed  out  to  four  decimal  place,  the  manufacturer  believed  the  program  accuracy  to  be 
somewhat  less  than  this. 

The  program  is  written  in  Pascal,  and  should  the  manufacturer  make  it  available, 
the  code  could  perhaps  be  modified  so  as  to  make  it  acceptable  to  conduct  simulations  of 
this  type.  Several  simulations  were  conducted  which  did  not  require  the  importation  of 
data  files.  These  runs  were  very  accurate  and  results  were  repeatable  during  multiple 
simulations.  Overall,  the  program  provides  very  accurate  orbital  positional  data  for 
predicting  and  visualizing  solar  system  events.  Without  modification,  however,  the 
program  is  not  sufficiently  precise  to  conduct  deflection  simulations  with  changes  of 
velocity  on  the  order  of  0.1  m/s. 

C.  DEFLECTION  RATIO  SENSITIVITY  TO  ORBITAL  PARAMETERS 

1.  Introduction 

Since  the  original  simulation  program  discussed  in  section  B,  above,  was  not 
satisfactory  to  test  the  deflection  hypothesis,  a  simulation  program  was  written  in 
MATLAB  which  permitted  the  application  of  a  given  change  in  velocity  to  a  variety  of 
hypothetical  asteroid  orbits.  The  MATLAB  scripts  used  are  contained  in  Appendix  C. 

The  change  in  velocity  was  applied  in  the  plane  of  the  asteroid’s  orbit.  The  direction  of 
the  change  in  velocity  was  applied  in  five  equal  increments  from  along  the  flight  path 
direction  to  the  anti-flight  path  direction.  As  the  answers  were  symmetrically  similar,  the 
directional  range  was  altered  to  five  equally  spaced  inputs  from  along  the  flight  path  to 
orthogonal  to  the  flight  path.  These  are  depicted  in  Figure  17. 
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Figure  17.  Deflection  Directions  With  Respect  to  the  Flight  Path  Vector. 

2.  Variation  in  Orbital  Period  for  a  Circular  Orbit 

The  first  sequence  of  tests  applied  a  0.135  m/s  change  in  velocity  to  an  asteroid 

with  a  circular  orbit  at  1  AU,  inclined  to  the  ecliptic.  The  Earth’s  orbit  was  modeled  as  a 
circular  orbit  in  the  ecliptic  at  1  AU.  The  simulation  calculated  the  orbital  elements  of 
the  asteroid’s  orbits  including  the  true  anomaly.  From  those  elements,  the  eccentric 
anomaly,  mean  anomaly,  and  mean  motion  were  calculated.  The  orbital  velocities  at  one, 
two,  three,  four,  five,  and  ten  years  prior  to  collision  were  calculated,  the  desired  change 
in  velocity  applied,  and  new  orbital  elements  calculated.  The  entire  procedure  was  then 
reversed.  Finally,  the  asteroid’s  position  at  the  original  time  was  determined,  compared 
to  the  Earth’s  position,  and  divided  by  the  Earth’s  radius  to  determine  a  miss  ratio  in 
Earth  radii.  As  expected,  the  asteroid  missed  by  two  Earth  radii  for  the  flight-path  and 
anti-flight  path  directions.  Results  show  that  the  miss  ratio  decreased  by  the  sine  of  the 
angle  between  the  flight  path  and  deflection  direction.  At  45  deg.,  the  miss  ratio 
decreased  to  1 .4  Earth  radii,  while  at  90  deg.,  there  was  no  deviation  detected.  The 
magnitude  of  the  deflection  distances  were  symmetric  about  90  deg. 
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3.  Constant  Radius  of  Perihelion 


To  demonstrate  the  accuracy  of  the  mathematical  solution  developed  in  Chapter 
IV ,  miss  ratios  were  calculated  for  an  asteroid  orbit  with  an  fixed  radius  of  perihelion  of  1 
AU.  Fixing  perihelion  distance  results  in  varying  eccentricity  and  semi-major  axis 
lengths.  Change  in  miss  ratio  with  change  in  eccentricity  and  change  in  semi-major  axis 
are  shov/n  in  Figure  18  and  Figure  19,  respectively.  Increasing  eccentricity  dictates  an 
in.ciease  in  semi-major  axis  for  a  constant  radius  of  perihelion.  As  such,  the  two  figures 
loc  i'  -i^eiy  much  the  same.  As  the  eccentricity  increased  from  a  minimum  value  of  0.0  to 
appi  oximately  0.2,  the  miss  ratio  increased  from  2  Earth  radii  to  approximately  4.2  Earth 
radii.  Over  this  same  variation,  the  semi-major  axis  increased  from  1  AU  to 
approximately  1 .25  AU.  As  shown  in  Eqn.  (24),  the  deflection  distance  varies  linearly 
with  the  asteroids  period.  But,  because  the  asteroids  period  varies  with  c^'^,  there  is  a 
relatively  large  increase  in  miss  ratio  or  increasing  values  of  eccentricity. 


i  ■4-4.5 
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I  *3-3 .5 
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Figure  18.  Miss  Ratio  in  Earth  Radii  vs.  Deflection  Angle  vs.  Eccentricity. 
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Figure  19.  Miss  Distance  in  Earth  Radii  vs.  Deflection  Angle  vs.  Semi-Major  Axis. 


4.  Variation  in  Eccentricity 

The  effect  of  varying  the  eccentricity  of  the  orbit  while  maintaining  a  constant 
semi-major  axis  length,  shown  in  Figure  20,  provided  several  interesting  results.  As  the 
eccentricity  increases  and  a  is  held  constant,  the  perihelion  distance  decreases.  The  line 
of  apsides  must  then  rotate,  increasing  the  true  anomaly,  to  keep  the  collision  distance  at 
1  AU.  The  true  anomaly  increases  approximately  5.7  deg.  for  each  0.1  increase  in 
eccentricity.  As  the  eccentricity  increases,  the  miss  distance  decreases.  The  primary 
cause  of  this  is  the  application  of  the  change  of  velocity  away  from  perihelion.  The 
greatest  deflections  are  achieved  by  applying  the  change  in  velocity  at  perihelion.  The 
least  efficient  deflections  are  at  aphelion.  In  this  case,  the  change  in  velocity  is  applied 
approximately  91  deg.  past  perihelion  for  eccentricity  of  0.01.  This  increases  to 
approximately  102  deg.  past  perihelion  at  e  =  0.2.  The  gradient  of  miss  ratio  with 
increase  in  eccentricity  increases  significantly  beyond  e  =  0.2.  This  shows  the  importance 
of  applying  a  given  change  in  velocity  at  perihelion. 
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Figure  20.  Miss  ratio  vs.  Deflection  Angle  vs.  Semi-Major  Axes. 

5.  Conclusions 

The  results  presented  above  are  consistent  with  the  mathematical  conclusions 
presented  in  Chapter  IV.  One  of  the  most  important  factors  in  maximizing  the  miss  ratio 
is  targeting  the  change  in  velocity  to  occur  at  perihelion.  This  is  especially  important  for 
high  eccentricity  orbit  with  a  semi-major  axis  close  to  that  of  the  Earth.  Much  larger 
changes  in  velocity  are  required  in  these  scenarios  the  further  from  perihelion  the  change 
is  to  occur. 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  SUMMARY 

The  issues  surrounding  EGA  mitigation  necessitate  a  wide  variety  of 
professionals,  each  with  differing  areas  of  expertise,  form  a  project  team  to  formulate  the 
most  feasible  solution.  This  paper,  for  the  Naval  Postgraduate  School  and  the  Space 
Warfare  Center,  is  the  genesis  of  another  small  portion  of  that  project  team.  Although  in 
its  infancy,  this  program,  because  of  the  unique  composition  of  its  members,  can 
contribute  unique  and  meaningful  contributions  to  that  effort. 

This  paper  serves  several  important  purposes.  First,  it  summarizes  a  few  of  the 
important  issues  involved  in  planetary  defense  and  serves  as  a  primer  to  helps  to  gain  an 
understanding  of  EGAs.  Secondly,  it  serves  to  develop  and  analyze  the  sensitivity  of 
deflection  angle  and  various  orbital  parameters  to  miss  ratio.  Most  importantly,  it 
identifies  areas  for  further  research  and  the  continuation  of  this  important  problem. 

Dance  of  the  Planets  is  not  sufficient  to  serve  as  a  model  for  the  purpose  of  testing 
deflection  hypotheses.  The  errors  induced  into  orbital  elements  or  vectors  when 
importing  data  are  of  sufficient  magnitude  to  prevent  accurate  analysis  of  results.  Should 
the  source  code  become  available,  increasing  the  precision  used  for  calculations  may 
solve  this  problem.  This  is  not  to  say,  however,  that  the  program  is  without  merit.  Used 
as  a  visualization  tool,  it  provides  an  extremely  flexible  and  accurate  model  with  which  to 
observe  intercept  geometry. 

The  variation  of  miss  ratio  with  deflection  angle,  eccentricity,  and  semi-major  axis 
changes  was  analyzed.  Deflection  along  the  flight  path  achieve  the  greatest  miss  ratios. 
Position  in  the  orbit  also  play  an  important  role  in  the  magnitude  of  the  miss  ratio.  For 
eccentric  orbits,  deflection  at  perihelion  is  important  as  miss  ratio  is  fairly  sensitive  to 
true  anomaly.  As  true  anomaly  at  deflection  increases,  the  miss  ratio  decreases.  This  a 
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direct  result  of  the  deflection  being  applied  at  further  angular  distances  from  perihelion 
resulting  in  less  of  a  change  in  period  for  each  asteroid  orbit. 

B.  FURTHER  RESEARCH  OPPORTUNITIES 

There  are  three  main  areas  in  the  asteroid  mitigation  problem  which  would  benefit 
from  additional  research.  First,  a  complete  engineering  and  mathematical  problem 
formulation  should  be  accomplished.  Expanding  the  scope  of  the  current  investigation  to 
include  a  review  of  the  state-of-the-art  in  change  in  velocity  delivery  options,  intercept 
trajectory  optimization,  and  weapons  technology. 

Second,  asteroid  orbit  determination  accuracy  versus  the  change  in  velocity 
required  to  mitigate  a  collision  should  be  researched.  The  accuracy  of  ephemeris  data  is 
extremely  important  when  calculating  required  changes  in  velocity  and  deflection  angle. 
The  relationship  between  the  two  should  be  determined  to  ensure  that  the  proper 
mitigation  steps  are  taken. 

Lastly,  the  formulation  of  the  sensitivity  model  developed  in  this  paper  should  be 
expanded  to  include  gravitational  perturbations  and  solve  the  problem  using  numerical 
integration  techniques.  A  high  fidelity  model  may  provide  additional  insights  into  the 
sensitivity  of  miss  ratio  to  deflection  angle. 


64 


LIST  OF  REFERENCES 


1 .  Morrison,  D.,  editor.  1 992.  The  Spaceguard  Survey:  Report  of  the  NASA 
International  Near -Earth-Object  Detection  Workshop.  (Pasadena:  Jet  Propulsion 
Laboratory). 

2.  Chapman,  C.  R.,  Harris,  A.  W.,  and  Binzel,  R.  1994.  Physical  properties  of  near-Earth 
asteroids:  Implications  for  the  hazard  issue.  In  Hazards  Due  to  Comets  and  Asteroids, 
ed.  T.  Geherels  (Tucson:  Univ.  of  Arizona  Press),  pp.  537-549. 

3.  Rahe,  J.,  Vanysek,  V.,  and  Weissman,  P.  R.  1994.  Properties  of  cometaiy  nuclei.  In 
Hazards  Due  to  Comets  and  Asteroids,  ed.  T.  Geherels  (Tucson:  Univ.  of  Arizona 
Press),  pp.  597-634. 

4.  Dance  of  the  Planets:  A  Dynamic  Model  of  the  Sky  and  Solar  System,  User’s  Manual. 
1995.  Applied  Research  &  Consulting,  Inc.  (Loveland:  Applied  Research  and 
Consulting). 

5.  Rabinowitz,  D.,  Bowell,  E.,  Shoemaker,  E.  M.,  and  Muinonen,  K.  1994.  The 
population  of  Earth-crossing  asteroids.  In  Hazards  Due  to  Comets  and  Asteroids,  ed. 

T.  Geherels  (Tucson:  Univ.  of  Arizona  Press),  pp.  285-312. 

6.  Shoemaker,  E.  M.,  Williams,  J.  G.,  Helin,  E.  F.,  and  Wolfe,  R.  F.  1979.  Earth¬ 
crossing  asteroids:  Orbital  classes,  collision  rates  with  Earth,  and  origin.  In.  Asteroids, 
ed.  T.  Gehrels  (Tucson:  Univ.  of  Arizona  Press),  pp.  253-282. 

7.  Shoemaker,  E.  M.,  Weissman,  P.  R.,  and  Shoemaker,  C.  S.  1994.  The  flux  of  periodic 
comets  near  the  Earth.  In  Hazards  Due  to  Comets  and  Asteroids,  ed.  T.  Geherels 
(Tucson:  Univ.  of  Arizona  Press),  pp.  313-335. 


65 


8.  Helin,  E.  F.  1982.  Earth-crossing  asteroids:  New  discoveries.  In  Sun  and  Planetary. 
System,  ed.  G.  Teleki  (D.  Reidel  Publishing  Company),  pp.  269-276. 

9.  Simonenko,  V.  A.,  Nogin,  V.  N.,  Petrov,  D.  N.,  Shubin,  O.  N.,  and  Solem,  J.  C.  1994. 
Defending  the  Earth  Against  Impacts  From  Comets  and  Asteroids.  In  Hazards  Due  to 
Comets  and  Asteroids,  ed.  T.  Geherels  (Tucson:  Univ.  of  Arizona  Press),  pp.  929- 
953. 

10.  Bate,  R.  R.  Mueller,  D.  D.,  and  White,  J.  E.  1971.  Fundamentals  of  Astrodynamics. 
(New  York:  Dover  Publications). 

11.  Wertz,  J.  R.  1978.  Editor,  Spacecraft  Attitude  Determination  and  Control.  Boston: 
Kluwer  Academic  Publishers. 

12.  Canavan,  G.  H.,  Solem,  J.  C.,  Rather,  J.  D.  C.,  eds.  1993.  Proceedings  of  the  Near- 
Earth-Object  Interception  Workshop.  (Los  Alamos:  Los  Alamos  National 
Laboratory). 

13.  Shoemaker,  E.  M.,  and  Helin,  E.  F.  1978.  Earth-approaching  asteroids:  Populations, 
origin  and  compositional  types.  NASA  Conf.  Public.  2053,  pp.  161-175. 


66 


APPENDIX.  MATLAB  SCRIPT  FILES 


This  appendix  contains  the  MATLAB  function  files  constructed  during  the 
research  of  the  asteroid  mitigation  problem.  Section  I  contains  the  scripts  which  calculate 
a  hypothetical  asteroid’s  position  and  velocity  vectors  based  on  the  current  position  and 
velocity  of  the  Earth.  Demo.m  calculates  all  data  by  accessing  individual  script  files.  In 
Demo.m  the  user  inputs  the  Earth’s  position  and  velocity  vectors,  rl  and  vl,  respectively. 
Position  is  in  1E6  km  and  velocity  is  in  km/sec.  For  an  elliptical  orbit,  the  user  inputs  into 
elip.m  the  desired  semimajor  axis  length,  eccentricity,  and  inclination.  For  a  parabolic 
orbit,  the  user  inputs  into  parab.m  the  desired  semi-latus  rectum,  and  inclination. 
Hyperbolic  orbit  script  files  were  not  developed.  When  demo.m  is  executed,  available 
outputs  include  the  asteroid’s  position  and  velocity  vectors  in  heliocentric  ecliptic 
coordinates  and  the  asteroid’s  orbital  elements.  Current  position  in  the  orbit  is  output  as 
true  anomaly. 

Section  2  calculates  the  miss  ratio,  the  ratio  of  the  difference  between  the  Earth’s 
position  and  the  asteroid’s  position  one  asteroid  orbit  following  the  imparted  change  in 
velocity  to  the  radius  of  the  Earth.  The  user  inputs  the  Earth’s  position  and  velocity 
vectors,  as  discussed  above,  into  Demol.m.  The  user  also  inputs  into  demo.m  the 
magnitude  of  the  change  in  velocity  to  be  applied,  normalized  for  one  orbit,  and  the 
number  of  orbits  prior  to  collision  the  change  in  velocity  is  to  be  applied.  In  elipl.m  the 
user  inputs  are  the  desired  semi-major  axis,  eccentricity,  inclination  for  asteroid  at  the  time 
of  Earth  impact.  The  program  then  calculates  the  remaining  orbital  parameters  for  the 
asteroid’s  orbit.  From  these,  the  eccentric  anomaly,  mean  anomaly,  and  mean  motion  are 
found  and  used  to  find  the  asteroid’s  position  at  some  specified  number  of  asteroid  orbits 
prior  to  collision.  The  change  in  velocity  is  then  applied  in  five  22.5  deg.  increments  from 
along  the  flight  path  to  90  deg.  to  the  flight  path.  New  values  for  mean  motion,  mean 
anomaly,  and  eccentric  anomaly  are  found  and  the  new  position  obtained  at  the  original 
time.  The  miss  ratio  is  then  calculated  to  determine  how  many  Earth  radii  the  new 
position  differs  from  the  old  position.  Outputs  include  new  velocity  vectors,  orbital 
elements,  true  anomaly,  and  miss  ratio  for  each  of  the  five  directions. 
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Summary  of  function  files: 


trunanom: 

True  anomaly  determination  from  r,v,e,n. 

incl; 

Inclination  determination  fi'om  h. 

Ian: 

Long,  of  Ascending  Node  determination  firom  n. 

aop: 

Argument  of  Perihelion  determination  fi'om  n,e,r. 

rot: 

Rotation  matrix  firom  from  perifocal  to  HCI.  Inputs  are 

Omega, i, omega. 

clip: 

Finds  r,v,nu,dcm  from  Earth’s  rmag,dcm,nu. 

parabola: 

Finds  r,v,nu,dcm  from  Earth’s  rmag,dcm,nu. 

earthorb: 

Calculates  Earth’s  r  vector  in  perifocal  and  HCI  from  a,ecc,dcm. 

cx,cy,cz: 

Compute  a  single  axis  direction  cosine  matrix. 

tlop; 

Compute  true  longitude  of  perihelion  from  e. 

deltav: 

Computes  the  nev/  velocity  vector  components  after  the  delta  v  is 
imparted  to  the  asteroid. 

elements: 

Returns  the  matrix  sets  of: 

el_prev:  orbital  elements  prior  to  deflection. 

el_next:  orbital  elements  after  deflection. 

rv_next:  position  and  velocity  vectors  at  original  impact 

time. 

elip  1 : 

Computes  the  asteroid’s  orbit  given  the  Earth’s  rmag,dcm,nu  and 
user  defined  values  for  a,e,i. 

elip2: 

Computes  the  asteroid’s  new  r,v  given  the  asteroid’s  position  at 
deflection,  r,  the  dcm  firom  the  orbital  elements,  the  position  in  the 
old  orbit,  ecc_a,a_a. 

function  1/2: 

Solve  Kepler’s  equation  for  E. 

param: 

Given  r,v  returns  a,ecc,Omega,incl,omega,nu,E,M,dcm,rmag. 

Variables: 

_ma^. 

r: 

Position  vector. 

v: 

Velocity  vector. 

e: 

Eccentricity  vector. 

n: 

Node  vector. 

h: 

Angular  momentum  vector. 

Omega: 

Longitude  of  the  ascending  node. 

i: 

Inclination. 

omega: 

Argument  of  perihelion. 

dcm: 

Direction  cosine  matrix  from  perifocal  to  heliocentric-echptic 
coordinate  system  for  Earth’s  orbit. 

nu: 

True  anomaly. 

a: 

Semi-major  axis. 

ecc: 

Eccentricity. 

dcm2: 

Direction  cosine  matrix  fi'om  perifocal  to  heliocentric-ecliptic 
coordinate  system  for  asteroid’s  orbit. 
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*_a 

*_old 

*_prev 

*_new 

*_next 

P 

r_P 

r_a 

E 

M 

n 


Earth  variable. 

Asteroid  variable. 

Asteroid  variable  at  time  of  impact. 

Asteroid  variable  at  some  time  prior  to  impact. 

Asteroid  variable  just  after  application  of  change  in  velocity. 
Asteroid  variable  in  new  orbit  but  at  time  of  original  impact. 
Semi-latus  rectum. 

Perihelion  distance. 

Aphelion  distance. 

Eccentric  anomaly. 

Mean  anomaly. 

Mean  motion. 
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L  ORBIT  DETERMINATION  AND  GRAPHING  SCRIPTS 


A.  CALCULATE  AND  PLOT  ORBITS:  FILE  ENTITLED  DEMO.M 


%  Given: 

%  First: 

%  Second: 
%  Third: 

%  a. 

%  b. 

%  c . 

%  Fourth: 
%  a. 

%  b. 


The  Earth's  position  and  velocity  vectors. 

Calculate  the  Earth's  orbital  elements. 

Plot  the  Earth's  orbit  wrt  the  mean  ecliptic. 

Choose  the  type  of  assailant  object  orbit: 

Circular  or  elliptical 

Parabolic 

Hyperbolic 

Choose  the  desired  values  for: 

Semi -major  axis  or  parameter 
Eccentricity 

New  inclination  with  respect  to  Earth's  inclination. 


clear 

%  Constants 

rtd=180/pi ; 
dtr=l/rtd; 

I  =  [10  0]; 

J  =  [0  10]; 

K  =  [0  0  1]; 

AU  =  1.4959965e8; 
TU  =  5.0226757e6; 
SU  =  29.784852; 
mu  =  1; 


%  Radians  to  degrees 
%  Degrees  to  radians 
%  Unit  vectors 


%  km 
%  sec 
%  km/sec 
%  km'^3/sec'^2 


%  Earth  values : 

rl=[-26  144.8  -.00038] ; 
vl=[-29.8  -5.376  -.000043]; 

r=rl*le6/AU; 
v=vl/SU; 

rmag=norm(r) ; 
vmag=norm{v) ; 

%  angular  momentum 
h=cross (r, v) ; 

hmag=norm(h) ; 

%  eccentricity  vector 

e=cross (v,h) -r/rmag; 
ecc=norm{e); 

if  ecc<le-10/  %  Check  if  eccentricity  is  zero 

e  =  [0  0  0]; 
ecc-0; 


%  Position 
%  Velocity 
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<i\o  o\<> 


end 


%  node  vector 
n  =  cross (K,h) ; 

nmag=norm{n)  ; 

if  nmag<le-10,  %  Check  if  inclination  is  zero 
n  =  [0  0  0]; 
nmag=0; 

end 

%  Find  the  true  anomaly  at  epoch 
nu  truanom (r , V,  e ,  n)  ; 

\  Compute  the  inclination 
inc  -■  incl(h) 

%  Compute  the  longitude  of  ascending  node 
Omega=lan (n) 

%  Compute  the  argument  of  perigee 
omega  =  aop(n,e,r) 

%  Compute  the  rotation  matrix 
dcm=rot (Omega, inc, omega) ; 

%  Determine  if  user  wants  another  orbit  calculated 
pit  =  input (’Do  you  want  to  calculate  a  new  orbit?  y/n  [y] :',’s'); 

if  isempty(plt) 
pit  =  ’y’; 

end 


pit  *y*, 

>•<*★**  **'*>**:(r-*-***'*-*-'55r*.********.***.5(fVlr****’******************'***'* 

: -loose  an  a,e  for  orbit  and  solve  for  nu__2  for  the  new  orbit 

f  priatf  ( '  \n‘Vca  may  now  calculate  a  new  intercept  orbit  typeAn'); 

typ:^  input  (  '■Do  you  want  to  plot  an  Ellipse,  Parabola,  or  Hyperbola?  e/p/h 
iz]  : ’ . ’s' ) ; 

if  isem.pty (type) 
type  =  *  e ’ ; 

end 


if  type  ==  ’e*,  %  Eliptical/Circular  Plot 

[r_new, v_new,nu_new,dcm2] =elip (rmag, dcm, nu)  ; 

elseif  type  ==  *p' ,  %  Parabolic  Orbit  Plot 

[r_new, v_new,nu_new] =parabola (rmag, dcm, nu) 

elseif  type  ==  *h’,  %  Hyperbolic  Orbit  Plot 

[r_new,  v__new,  nu_new]  =hyperb  ( rmag ,  dcm ,  nu ) 
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end  %  end  plotting  of  second  orbit 
r_old  =  rl 

r_new  =  r_new  *  AU/le6 

v_old  =  vl 

v_new  =  v_new*SU 

del_v  =  abs (norm (v_old~v_new) ) 

end 


Calculate  the  Earth's  orbit 


p  =  hmag'^2/mu; 
a  =  p/  (l“ecc'^2)  ; 
rp  =  p/ (l+ecc) ; 
ra  =  p/ (1-ecc) ; 


semilatus  rectum 
semimajor  axis 
perigee 
appogee 


[x,y, z , XX, yy, zz] =earthorb (a, ecc, dcm) ; 


%  Plot  the  Earth's  orbit  with  respect  to  the  ecliptic 
hold  on; 

plots {xx,yy, zz, 'm' ) ;  %  orbit 


xyz_r=r;  %  position 

plots  (xyz_r  (1)  ,xyz_r  (2)  ,xyz_r  (S)  ,  'g+'  )  ; 

if  type  ==  'e',  %  Eliptical/Circular  Plot 

legend ( ' b ' , ' Asteroid ' ' s 

orbit '  ,  ' or’ , ' Perihelion' , ' *r ' , 'Aphelion' , ' +g’ , ' Current  Position' , . . . 
'm' , ' Earth' ' s  orbit ' ) 


elseif  type  ==  'p'  |  'h',  %  Parabolic  Orbit  Plot 

legend ( ' b ' , 'Asteroid* ' s  orbit  * , *  or  * , ' Perihelion' , ' +g' , ' Current 
Position' , . . . 

'm' ,' Earth' * s  orbit’) 

end 


xyz_per=dcm*  [rp; 0; 0]  ;  %  perihelion 

plots  (xyz_per  (1)  ,xyz__per  (2)  ,xyzj)er  (S)  ,'mo'); 

xyz__apo=dcm*  [-ra;  0;  0]  ;  %  aphelion 

plots (xyz_apo (1) ,xyz_apo (2) ,xyz_apo (S) , 'm* ' )  ; 

%  end  plotting  of  second  orbit 

%  Axis  plotting  statements. 

1=2; 

ds=.2; 

A=  [0  0  0;1  0  0]  ; 

B=[0  0  0;0  1  0]  ; 

C=[0  0  0;0  0  1]  ; 
plots  (A,B,C,  '  w'  )  ; 

text (1+ds, 0, 0, 'X' , ' horizontalalignment ' , 'center' )  ; 
text ( 0 , 1+ds , 0 , ' Y' , ' horizontalalignment  * , ' center'  )  ; 
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text (0,0, 1+ds , ' Z ' , ’ horizontal alignment ' , ’ center' ) 

view{viewmtx(135,30,25) ) ; 

title (' General  View'); 

axis (' square '); axis  off; 

axis ([-a,a,-a,a,-a,a]); 

end; 

View  radio  buttons 
txt_view=uicontrol (gcf , . . . 

' Style ' , ' text ’ , . . . 

'String' , 'Asteroid  Orbit  View' , , . . 

'Units' , 'normalized' , . . . 

'Position' ,  t.05, .87, .20, .03] ) ; 
view_gen=uicontrol (gcf, . . . 

' Style' , 'radio' , . . . 

'String' , 'Standard' , . . . 

'Units' , 'normalized' , . . . 

' Position' ,  [ . 05,  .  80 ,  .20, . 05] , . .  . 

'Value' , 1, . . . 

' CallBack ' ,  [ . . . 

' title (* 'General  View' 


set (view_gen. 

' 'Value 

'  M)  , 

set (view_x, * ' 

Value*  * 

,0)  ,  '  . 

set (view_y, ' ' 

Value '  ' 

,0)  ,  '  . 

set {view_z, ' * 

Value '  ' 

,0)  ,  ’  . 

set  (view_n,  '  * 

Value*  * 

,0)  ,  '  . 

'view (135, 30) ' ] ) ; 
view_x=uicontrol (gcf, . . . 

' Style ' , ' radio* , . . . 

' String ' , ' Down  X  Axis ' , . . , 

'Units ' , 'normalized* , . . . 

'Position' ,  [.05,  .75, .20, .05] , . . , 

' CallBack' , [ . . . 

' title ( ' 'View  Looking  Down  +X  axis '*),'... 

'  set  (view__gen,  '  'Value '  '  ,  0)  ,  *  .  .  . 

' set (view_x, ’ 'Value 
' set {view_y, * 'Value ’  ' , 0) , '  .  .  . 

' set ( view_z , ' ' Value ' ' , 0 ) , ' . . . 

' set (view_n, ' 'Value ' * , 0) , * . . . 

'view( [1,0,0] ) '] ) ; 

view_y=uicontrol (gcf,  .  .  . 

' Style ' , ' radio ' , . . . 

' String ' , ' Down  Y  Axis ' , . . . 

'Units ' , 'normalized' , . . . 

' Position' ,  [ . 05,  .70, .20, . 05] ,  .  .  . 

' CallBack' ,  [ . . . 

' title (' 'View  Looking  Down  +Y  axis*'), 
*  set (view_gen, ' 'Value ' ' , 0) , ' . . . 

'  set  (view__x,  '  'Value  *  ’  ,  0)  ,  '  . .  . 

'  set  (view__y,  '  'Value '',!),'... 

' set (view_z , ' 'Value ' ' , 0) , ' . . . 

' set {view_n, * 'Value'  ' , 0) , ' .  .  . 

' view(  [0,1,0])']); 
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view_z=ui control (gcf , , . . 

'Style' , 'radio' , . . . 

'String* , 'Down  Z  Axis' , . . . 

'Units' , 'normalized' , . . . 

' Position' ,[.05,  .65, .20,  .05],... 

'CallBack' ,  [. . . 

' title  ( '  'View  Looking  Down  +Z  axis ''),*... 
' set (view_gen, ’ 'Value ’ '  ,  0)  ,  '  . .  . 


set  (view_ 

'Value' 

'  ,0) 

set  (view_ 

-Y/  ' 

'Value' 

'  ,0) 

set  (view__ 

J-'  ' 

' Value ' 

M) 

set (view 

n,  ' 

'Value' 

'  ,0) 

'view( [0,0,1]  )  ']  )  ; 

xyz_norm=dcm2  * [ 0 ; 0  ;  1  ]  ; 
view_n=ui control (gcf , . . . 

' Style ' , *  radio  * , . . . 

'String' , 'Down  Orbit  Normal' , . . . 

'Units' , 'normalized' , . . . 

' Position' ,[.05, .60, .20, .05],... 

'CallBack' ,  [. . . 

' title (' 'View  Looking  Down  Orbit  Normal''),'... 
'  set  (view__gen,  '  'Value '  '  ,  0)  ,  '  .  .  . 


set  (view_x,  ' 

'Value' 

'  ,0) 

set (view_y, ' 

■Value' 

’  ,0) 

set (view_z, ' 

'Value' 

\0) 

set (view  n,  ' 

'Value' 

M) 

'view(  [xyz_norm(l)  ,xyz_norm(2)  ,xyz_norm(3)  ])*]); 

%  Zoom  Slider 
zoom=uicontrol (gcf,  .  .  . 

' Style' , ' slider' , , . . 

'Units ' , 'normalized' , . . . 

' Position' ,  [ .  05,  ,  15, .20, . 05]  ,  .  . . 

'Min’ , .5, 'Max' ,5.5, 'Value' ,1, . . . 

'CallBack' , [ . . . 

' set (zoom_cur, ' ' String' ' ,num2str (get (zoom, '  'Val*  *))),’,. 
' set (gca, ' 'Xlim' * , [- 

get (zoom, ' 'Value' ' ) ,get (zoom, * 'Value' ')],',... 

' ' ' Ylim' ' ,  [-get (zoom, ' 'Value ' * ) ,get (zoom, ' 'Value*  *)],',. 
' ' ' Zlim'  ' ,  [-get (zoom, ' 'Value ' ' ) ,get (zoom, ' 'Value*  * )  ]  )  '  ]  ) 
zoom_label=uicontrol (gcf, . . . 

' Style’ , ' text ' , . . . 

'Units ' , 'normalized' , . . . 

' Position' ,  [ , 05,  .205 ,  .20 , . 03] ,  . . . 

' String* , ' Zoom’ ) ; 
zoom_cur=uicontrol (gcf, . . . 

'Style' , 'text' , . . . 

'Units' , 'normalized' , . . . 

*  Position* ,  [.095, .115,  .11,  .03]  ,  ... 

' String' , num2str (get (zoom, 'Value ' ) ) ) ; 
zoom_in=ui control (gcf, . . . 

' Style ' , ' text ' , . . . 

'Units ' , 'normalized' , . . . 


o\o  o\o 


'Position' ,  [.05, .115, .04, .03] , . . . 
' String' , ' In' ) ; 
zoom_out=uicontrol (gcf , . . . 

' Style ' , ' text ' , . . . 

' Units ' , ' normalized ' , . . . 
'Position' ,  [.21, .115, .04, .03] , . . . 
'String' , 'Out' ) ; 

%  Print  push  button 
prt 


B.  AOP 

% 

* 

% 

Calculation  of  the  argument  of 

* 

% 

perigee,  in  degrees 

* 

* 

% 

LCDR  Wade  Knudson 

★ 

Naval  Postgraduate  School 

* 

o, 

*0 

Orbit  Visualization  Routines 

★ 

% 

★ 

% 

Inputs:  n,e,r 

★ 

% 

* 

% 

Output :  omega 

* 

% 

* 

% 

Files  called:  tlop 

* 

“o 

* 

% 

[omega]  =  aop{n,e,r) 

* 

% 

* 

%- 

C.  ARGUMENT  OF  PERfflELION 

function  [omega]  =  aop(n,e,r) 

nmag=norm(n) ; 
ecc=norm(e) ; 
rmag=norm(r)  ,- 
I  =  [10  0]; 

J  =  [0  10]; 

K  =  [0  0  1]; 
rtd=180/pi; 


Check  to  see  if  the  orbit  is  circular  (e=0)  and  in  the  ecliptic  (n=0) 
if  so,  then  omega  is  not  used;  use  true  longitude  at  epoch  (Vallado  p 
139) 


if  nmag  ==  0  &  ecc  ==  0, 
omega  =  0; 
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%  Check  to  see  if  the  orbit  is  elliptical  (e>0)  and  in  the  ecliptic  (n=  0)  ; 
%  if  so,  then  omega  is  not  used;  use  true  longitude  of  periapsis  (Vallado 
p,  136) 

elseif  nmag  ==  0  &  ecc 

om__true  =  tlop(e); 

%fprintf ( *  This  is  an  elliptical  orbit  in  the  ecliptic  An’) 

%fprintf ( ' The  inclination  is  zero  An’) 

%fprintf  (  An’ ) 

%fprintf ( ’ The  angle  between  Aries  and  periapsis  is  called  \n*) 
%fprintf ( ' the  true  longitude  of  periapsis,  and  is  %3.1f 
deg  An'  ,om_true) 

omega  =  om_true; 


%  Check  to  see  if  the  orbit  is  a  circular  (e=0)  and  inclined  to  the 
ecliptic 

elseif  nmag  ~=  0  &  ecc  ==  0, 
if  dot (r, K) <0, 

omega  =  360  -  acos (dot (n, r) / (nmag*rmag) ) *rtd; 

else 

omega  =  acos (dot (n, r) / (nmag*rmag) ) *rtd; 

end 


%  All  remaining  cases 
else 

if  dot(e,K)  <  0, 

omega  =  360  -  acos (dot (n, e) / (nmag*ecc) ) *rtd; 
else 

omega  =  acos (dot (n, e) / (nmag*ecc) ) *rtd; 

end 


end 

D.  CX 

function  f=cx (angle) 

f=[l  0  0;  0  cos (angle)  sin (angle) ;  0  -sin (angle)  cos (angle) I ; 

KCY 

function  f=cy (angle) 

f= [cos (angle)  0  -sin (angle) ;  010;  sin (angle)  0  cos (angle)]; 

R  CZ 

function  f=cz (angle) 

f= [cos (angle)  sin(angle)  0;  -sin(angle)  cos (angle)  0;  001]; 
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G.  EARTH  ORB 


Calculation  of  the  Earth's  orbit 

LCDR  Wade  Knudson 

Naval  Postgraduate  School 

Orbit  Visualization  Routines 

Inputs:  a,ecc,dcm 

Output:  x,y, z:  perifocal  coordinates 

xx,yy,zz:  heliocentric  coordinates 

Files  called:  None 

[x,y, z,xx,yy, zz]  =  earthorb (a, ecc^dcm) 


★ 

★ 

* 

★ 

★ 

* 

* 

★ 

★ 

* 

* 


function  [x,y, z,xx,yy, zz]  =  earthorb (a, ecc , dcm) 


m=400; 

E=linspace(0,2*pi,m) ; 

x=-a*ecc+a*cos (E) ;y=a*sqrt (l-ecc^2) *sin(E) ;  %Battin  p.  158 
z=zeros (l,m) ; 
for  ]c=l:Tn 

point =dcm* [x{k) ;y (k) ;z (k) ] ; 

XX (k) =point (1) ;yy (k) =point (2) ; zz (k) =point (3) ; 

end; 


H.  FJ  IPTICAL  ORBIT  CALCULATIONS 


Calculation  of  the  Eliptical/circular  orbit 

LCDR  Wade  Knudson 

Naval  Postgraduate  School 

Orbit  Visualization  Routines 

Input  s :  rmag , dcm , nu 

User  inputs:  a,ecc,incl 

Output:  r,v,  and  nu  for  the  new  orbit 

Files  called:  cx,c2 


* 

* 

* 

* 

★ 

* 

* 

* 

* 

•k 

* 

* 

★ 

* 

★ 
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o\o  o\o  o\o 


[r,v,nu_new]  =  elip (rmag, dcm, nu) 


* 


* 


function  [r, v,nu_new,dcm2]  =  elip (rmag,dcm,nu) 

%  Constants 

rtd  =  180/pi; 
dtr  =  l/rtd; 
mu  =  1; 

%a  =  input ('Type  in  the  semi-major  axis') 

%  If  ecc  ==  0,  then  'a'  must  equal  rmag. 

%  If  ecc  0,  then  rp  <  a  <  ra. 
a=2.032; 
ecc=.651; 
newinc=5 . 46  ; 

if  ecc  ==  0,  %  If  circular  then  no 

%  rotation  to  argument 
nu_new  =0;  %  of  periapsis  (nu_new) 

a  =  rmag; 
p  =  a; 

else 

p  =a* (l-ecc^2) ; 

nu__new=acos  (  (p/rmag-l) /ecc)  *rtd;  %  true  anomaly 

end 

rp  =  p/ (1+ecc) ; 
ra  =  p/ (1-ecc) ; 

%  Perifocal  coordinates 

r  =  [rmag*cos (nu_new*dtr)  rmag*sin {nu_new*dtr)  0]  ; 

V  =  sqrt (mu/p) * [-sin (nu_new*dtr)  ecc+cos (nu_new*dtr)  0]  ; 


%  periapsis 
%  apoapsis 


%  rotate  the  new  orbit  to  the  line  of  node  (-nu) ,  then  to  desired 
inclination, 

%  and  finally  rotate  nu_2  to  get  epoch  to  r. 

%  Asteroid' s  dcm  matrix 

dcm2=  dcm*cz (-nu*dtr) *cx (-newinc*dtr) *cz (nu_new*dtr) ; 
m=200; 

E=linspace ( 0 , 2*pi ,m) ; 
x=-a*ecc+a*cos (E) ; 

y=a*sqrt  (l-ecc'^2)  *sin(E)  ;  %  Battin  p.  158 

z=zeros (l,m) ; 
for  k=l:m 

point =dcm* [x(k) ;y (k) ;z (k) ] ; 

XX (k) =point (1) ;yy (k) =point (2) ;zz (k) =point (3) ; 
point2=dcm2* [x (k) ;y (k) ; z (k) ] ; 

XXX (k) =point2 (1) ;yyy (k) =point2  (2) ;zzz (k) =point2 (3) ; 


79 


end 


%  Plot  the  asteroid' s  new  orbit 
hold  on; 

plots (xxx,yyy,zzz, 'b' ) ;  %  orbit 

xyz_per=dcm2* [rp;0;0] ;  %  perihelion 

plots (xyz_per (1) ,xyz_per(2) ,xyz_per{S) , 'bo' ) ; 

xyz_apo=dcm2* [-ra;0;0] ;  %  aphelion 

plots (xyz_apo ( 1 ) , xyz_apo (2) , xyz_apo ( S ) , ' b* ' ) ; 

xyz_r= (dcm2*r' ) ' ;  %  position 

%  rvector 

plots (xyz_r (1) ,xyz_r (2) ,xyz_r (S)  ,'g+'); 

legend ( 'Asteroid' 's  orbit' , 'Perihelion' , 'Aphelion' ) 

%  Calculate  the  new  r,v  vectors 
r  =  xyz_r  ; 

V  =  (dcm2*v' ) ' ; 


I. 

%-- 

INCLINATION 

% 

Q, 

•o 

% 

Calculation  of  the  inclination  in  degrees 

★ 

* 

% 

* 

% 

LCDR  Wade  Knudson 

* 

o, 

"o 

Naval  Postgraduate  School 

* 

% 

Orbit  Visualization  Routines 

* 

% 

★ 

o, 

“o 

Inputs :  h 

* 

% 

*• 

o, 

'0 

Output :  inc 

★ 

% 

* 

Q, 

*0 

Files  called:  none 

* 

'o 

% 

[inc]  =  inc(h) 

* 

* 

% 

★ 

%--- 

function  inc  =  inc(h) 

K-[0  0  1]  ; 
rtd=180/pi; 
hmag=norm (h) ; 

inc  =  acos (dot (h, K) /hmag) *rtd; 
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J.  LONGITUDE  OF  ASCENDING  NODE 


Calculation  of  the  longitude  of 
the  ascending  node,  in  degrees 


★ 

* 

* 

★ 


LCDR  Wade  Knuds on  * 

Naval  Postgraduate  School  * 

Orbit  Visualization  Routines  * 

* 

Inputs:  n  * 

* 

Output :  Omega  * 

* 


Files  called:  none 
[Omega]  =  lan(n) 


function  [Omega]  =  lan(n) 

nmag=norm(n) ; 

I  -  [10  0]; 

J  =  [0  10]; 
rtd=180/pi; 

if  nmag  ==  0,  %  Provides  check  for  inclination 

Omega  =0;  %  If  no  inclination  then  no  Omega 

else  %  If  some  inclination  then 

%  find  Omega 

if  dot (n, J)  <0, 

Omega  =  360  -  acos (dot (n, I) /nmag) *rtd; 
else 

Omega  =  acos (dot (n, I) /nmag) *rtd; 

end 

end 


K.  PARABOLA 


3-D  plot  of  the  parabolic  orbit  which  intercepts 
the  Earth’s  position  at  a  given  time 

LCDR  Wade  Knudson 

Naval  Postgraduate  School 

Orbit  Visualization  Routines 

Inputs:  rmag,dcm,nu 

User  input:  p,incl 

Output:  r,v,  and  nu  for  the  new  orbit 


* 

* 

* 

* 

* 

* 

* 

★ 

* 

★ 

* 

* 
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Files  called:  cx,cz 


[r,v,nu_new]  =  parabola (rmag,dcm,nu) 


* 

* 

* 

* 


* 


function  [r,v,nu_new]  =  parabola (rmag, dcm, nu) 

%  Constants 

rtd  =  180/pi; 
dtr  =  1/rtd; 
mu  =  1; 

P=1.9; 

newinc=90; 

ecc=l; 


nu_new=acos ( (p/rmag-l) /ecc) *rtd;  %  true  anomaly 

rp  =  p/{l+ecc)  %  perihelion 

%  Perifocal  Coordinates 

r  =  [rmag*cos {nu_new*dtr)  rmag*sin {nu_new*dtr)  0]  ; 

V  =  sqrt (mu/p) * [-sin (nu_new*dtr)  ecc+cos (nu_new*dtr)  0] ; 


%  rotate  the  new  orbit  to  the  line  of  node  (-nu),  then  to  desired 
inclination, 

%  and  finally  rotate  nu_2  to  get  epoch  to  r. 

dcm2=  dcm*cz (-nu*dtr) *cx ( -newinc*dtr) *cz (nu_new*dtr) ; 


m=200; 

f=l inspace (-pi+.l,pi- .l,m) ; 
xl  =  (p*cos (f ) ) ;yl  =  (p*sin(f)); 
x2  =  (1+cos (f ) ) ;y2  =  (l+cos(f)); 

X  =  xl./x2;y  =  yl./y2;  %  Battin  p.  158 

z=zeros (l,m) ; 

for  k=l:m 

point =dcm* [x (k) ;y (k) ; z (k) ] ; 

XX (k) =point (1) ;yy (k) =point (2) ; zz (k) =point (3 ) ; 
point2=dcm2* [x(k) ;  y{k) ;  z (k) ] ; 

xxx(k) =point2 (1) ;yyy(k) =point2 (2) ;zzz (k)=point2 (3) ; 

end 

%  Plot  the  asteroid's  new  orbit 
hold  on; 
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plot 3 (xxx,yyy, zzz, ’b' ) ; 


orbit 


xyz_per=dcm2* [rp; 0 ; 0] ;  %  perihelion 

plots  (xyzjper  (1)  ,xyz_per  (2)  ,xyz_per  (3)  ,*bo'); 

xyz_r= (dcm2*r  ’  )  '  ;  %  position 

%  rvector 

plots  (xyz_r  (1)  ,xyz_r(2)  ,xyz_r  (3)  ,  ’g+* )  ; 

%  Calculate  the  new  r,v  vectors  and  orbital  elements, 
r  =  xyz_r 
V  =  (dcm2*v' ) ' 

L.  PRINT 

%  Puts  a  print  push  button  on  a  figure, 
function  prnbut  =  prt 

prt=uicontrol (gcf , . . . 

' Style ' , 'push' , . . . 

'Units ' , 'normalized' , . . . 

'Position' ,  [.05, .02,  .2  0, .05]  , . .  . 

' String ' , ' Print ' , . . . 

' CallBack' , 'print ' ) ; 


M.  ROTATION  DCM 


Calculation  of  the  rotation  matrices  for  a  3-1-3 
rotation  to  convert  from  perifocal  coordinates  to 
heliocentric -ecliptic 

LCDR  Wade  Knudson 

Naval  Postgraduate  School 

Orbit  Visualization  Routines 

Inputs :  Omega , inc , omega 

Output :  dcm 

Files  called:  cx,cz 

[dcm]  =  rotmat (Omega, inc, omega) 


* 

* 

★ 


* 

★ 


★ 

* 

* 

* 

★ 

★ 

* 

* 

* 


* 

* 


* 


function  [dcm]  =  rot (Omega, inc, omega) 


dtr=pi/l80; 

0=0mega*dtr; 
i=inc*dtr ; 


83 


o=omega*dtr; 


dcm  =  cz (-0) *cx (-i) *cz (-o) ; 


N.  TRUE  LONGITUDE  OF  PERIHELION 


True  longitude  of  periapsis 

Used  for  elliptical  orbits  in  the  reference  plane 

LCDR  Wade  Knudson 

Naval  Postgraduate  School 

Orbit  Visualization  Routines 

Inputs:  e 

Output:  True  longitude  of  periapsis 
Files  called:  none 
[om_true]  =  tlop{e) 


★ 

★ 

★ 

* 

★ 

•k 

k 

k 

k 

k 

k 

k 

k 

k 

k 


function  om_true  =  tlop(e) 

rtd  =  180/pi; 

I  =  [10  0]; 

J  =  [0  10]; 
etnag  =  norm  ( e )  ; 

if  dot(J,e)  <  0, 

om_true  =  360  -  acos (dot (I , e) /emag) *rtd; 

else 

om_true  =  acos (dot (I, e) /emag) *rtd; 

end 


O.TRUE  ANOMALY 


Calculation  of  the  true  anomaly  at  epoch, 

LCDR  Wade  Knudson 

Naval  Postgraduate  School 

Orbit  Visualization  Routines 

Inputs:  r,v,e 

Output :  nu 

Files  called:  none 


in  degrees 


* 

★ 

* 

* 

* 

* 

* 

k 

k 

k 

k 

k 

k 
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[nu]  =  truanom(r, V, e,n) 


function  [nu]  =  truanom (r , v, e , n) 

rmag=norm (r) ; 
nmag=norm (n) ; 
ecc=norm(e) ; 
rmag=norm (r) ; 

I  =  [10  0]; 

J  =  [0  10]; 

K  =  [0  0  1]; 
rtd=180/pi; 

%  First  check  for  circular  orbit  in  the  ecliptic.  If  not,  skip, 
if  ecc  ==  0  &  n  ==  0, 
if  dot(J,r)<0, 

lamda__t  =  360  -  acos  (dot  ( I ,  r) /rmag)  *rtd; 
else 

lamda_t  =  acos (dot (I , r) /rmag) *rtd; 
end 

nu  =  lamda_t; 

%fprintf ( *  This  is  a  circular  orbit  in  the  ecliptic.  \n*) 

%fprintf ( *  The  longitude  of  ascending  node  is  undefined. \n' ) 

%fprintf ( ' The  inclination  is  zero,\n') 
fprintf ( ' \n' ) 

%fprintf ( *  The  angle  between  aries  and  the  position  of  the  body\n») 
%fprintf ( ' is  called  the  true  longitude  at  epoch,  and  is  %3.1f  deg,\n*,nu) 


%  Next  check  for  circular  inclined  orbit.  If  not,  skip, 
elseif  ecc  ==  0  &  nmag  0, 

if  dot(r,K)<0, 

u  =  36  0  -  acos  (dot  (n,  r)  /  (ntnag^rmag)  )  *rtd; 
else 

u  =  acos (dot (n,r) / (nraag*rmag) ) *rtd; 
end 

nu  =  u; 

fprintf (' This  is  a  circular,  inclined  orbit.  \n’) 
fprintf ( * \n' ) 

fprintf (* The  argument  of  perigee  is  known  as  the  \n’) 
fprintf (* argument  of  latitude  at  epoch,  and  is  %3.1f  deg.\n',u) 

else 

if  dot(r,v)  <  0, 

[nu]  =  360  -  acos (dot (e, r) / (ecc*rmag) ) *rtd; 
else 

[nu]  =  acos (dot (e , r) / (ecc*rmag) ) *rtd; 
end 

%fprintf ( ’ This  is  an  elliptical  orbit.  \n') 
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o\o  o\o 


if  nmag  ==  0  &  ecc  ~=  0, 

%fprintf ( ’The  orbit  is  in  the  ecliptic . \n' ) 

else 

%fprintf ( ' The  orbit  is  inclined  to  the  ecliptic . \n’ ) 

end 

fprintf ( ' \n' ) 

fprintf(’The  true  anomaly  at  epoch  is  %3.1f  deg.\n’,nu) 
end 
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n.  SENSITIVITY  SCRIPTS 


A.  MISS  RATIO  CALCULATIONS:  DEMOl 

clear  all 

rtd=180/pi; 

dtr=l/rtd; 

I  =  [10  0]; 

J  =  [0  10]; 

K  =  [0  0  1]; 

AU  =  1.4959965e8;  %  km 

TU  =  5.0226757e6;  %  sec 

SU  =  29.784852;  %  km/sec 

mu  =  1;  %  km''3/sec'^2 

%r=  [-26  144.8  - . 00038]  ; 
lrv=[-29.8  -5.376  -0.000043]; 


r=[AU  0  0]/le6; 
v=[0  SU  0]  ; 
num_orbs  =  1 ; 


%  Calculate  the  orbital  elements  for  the  Earth  from  the  r/v  vectors; 

%  Input  are  the  position  and  velocity  vectors 
%  Outputs  are  shown: 

[a_e , ecc_e , 0_e , inc_e , o_e , nu_e , E_e , M_e , dcm_e , rmag_e]  =  param ( r , v) ; 

S'  S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S-S'S'S'S'S'S'S' 
■o  'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'S'o'o'o'o'o'o'o'o'o'o'o'o'o'o'e'o'o'o'o'o'o'o'o'o'o'b'o'o'o'o'o'o'o'o'o'S'b'o'o'o'o'o'o'o'o'o 

% 

%  choose  an  a,e  for  orbit  and  solve  for  nu_2  for  the  new  orbit 


%  Determine  if  user  wants  another  orbit  calculated 
%plt  =  input ('Do  you  want  to  calculate  a  new  orbit?  y/n  [y]  :',*s'); 

%if  isempty(plt) 
pit  =  'y' ; 

%end 


if  pit  ==  'y' , 

%fprintf  ( '  \nYou  may  now  calculate  a  new  intercept  orbit  typeAn'); 

%type  =  input ('Do  you  want  to  plot  an  Ellipse,  Parabola,  or  Hyperbola? 
e/p/h  [e] : ' , ’ s ' ) ; 

%if  isempty (type) 
type  =  ' e ' ; 

%end 
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o\o  oV>  o\o 


Calculate  new  r/v  vectors  based  on  inputting  desired  values  of 
a,e,inc.  Omega, omega,  and  nu  are  not  variables,  but  fixed,  based 
on  the  Earth's  orbit. 

if  type  ==  »e',  %  Eliptical/Circular  Plot 

[r_new,  v_new,nu_new,dcm2,ecc]  =elipl  {rmag__e, dcm_e,nu_e)  ; 

elseif  type  ==  'p',  %  Parabolic  Orbit  Plot 

[r_new,  v__new,nu_new]  =parabola  (rmag,dcm,nu) 

elseif  type  ==:  'h',  %  Hyperbolic  Orbit  Plot 

[r_new,  v_new,  nu_new]  =hyperb  (rmag,  dcm,nu) 

end  %  end  plotting  of  second  orbit 

%  Compute  the  orbital  parameters  of  the  new  orbit: 
global  compare 

[ a_a ,  ecc_a ,  0_a ,  inc_a ,  o_a ,  nu_a ,  E_a ,  M_a ,  dcm_a ,  rmag_a]  =  param  ( r__new ,  v_new)  ; 

%  Compute  the  period  of  the  asteroid's  orbit 
mu_s  =  1 ; 

TP  =  2*pi*(a__a)"l.5/sqrt{mu_s)*TU/3600/24/365.25;  %  In  Earth  years 

%  Input  here  how  many  asteroid  orbits  to  go  through.  You  should  go  through 
%  at  least  one  asteroid  orbital  period  to  ensure  sufficient  time  for 
deflection, 

n_a=  sqrt  (mu_s/a_a^3)  ;  %  mean  motion 

t=-TP*num_orbs*365 .25*24*3600/TU;  %  convert  years  to  canonical 

global  t 

%  Compute  the  required  delta  V  required  along  the  flight  path 
del_tot=0 . 000135/num_orbs ; 

%  Compute  mean  anomaly  for  some  past  time. 

M_prev  =  M_a  +  n_a*t; 

M_prevl  =  M_prev  -  f ix (M_prev/2/pi) *2*pi  +  2*pi; 
global  M_prev  ecc_a 


%  Previous  eccentric  anomaly 
E_prev  =  f zero (' fun' ,M_prev) ; 

%  Compute  previous  true  anomaly 
nu_prev  =  acos ( {ecc_a-cos (E_prev) ) / (ecc_a*cos (E_prev) -1) ) *rtd; 

%  Compute  rmag  for  previous  time 
rmag__prev  =  a_a*  (l-ecc_a*cos  (Ej)rev}  )  ; 

%  Compute  the  r/v  at  the  previous  time 


[r_prev  v_j)rev]  =elip2  (rmag_prev,  dcm_a,  nu_prev,  ecc_a,  a_a)  ; 

%  [a_a_p ,  ecc_a_p ,  0_a_p ,  inc_a_p ,  o_a_p ,  nu_a_p ,  E_a_p,M_a_p ,  dcm_a_p]  = 
param{r_j}rev,  v_prev)  ; 


%  Compute  the  new  matrix  of  delta_v’s 

[el_prev  el_next  rv_next]=deltav(r_prev,v_prev,del_tot,nu_new); 

delta_rv= zeros (6,6)  ; 
for  s=l:6 

delta_rv(s,  1 :3)  =  (r__e  -  rv_next(s,l:3))*le6; 
delta_rv(s,4 ;6)  =  y_prev  -  rv_next (s , 4 :6)  ; 

mag_del_rv  (s ,  :)  =  [norm  (delta__rv  (s ,  1 : 3)  )  norm(delta_rv(s,  4  :  6)  )  ] 

end 

time_ratio=  [nu_new  a_a  ecc_a  (mag_del_rv  (2  : 6 , 1) /6378 . 145)  *  ]  ; 
fprintf(’\n  nu  a  e  0.0  22,5  45  67.5 

90\n») ; 
fprintf ( ' \n 

%7,4f%7.4f%7.4f%7.4f%7.4f%7.4f%7.4f%7.4f%7.4f\n\n* , [time_ratio] ) ; 

fprintf ( ’ \n' ) ; 
keyboard 

B.  DELTA  V  DETERMINATION 

function  [el__j>rev,  el_next ,  rv_next]  =  deltav  (rl ,  vl ,  dv_tot ,  nu__new) 

%  Constants 
rtd=180/pi; 
dtr=l/rtd; 

I  =  [10  0]; 

J  =  [0  10]; 

K  =  [0  0  1]; 

AU  =  1.4959965e8; 

TU  =  5.0226757e6; 

SU  =  29.784852; 
mu  =  1; 

%  Orbit  values ; 
r=rl*le6/AU; 
v=vl/SU; 

rmag=norm(r) ; 
vmag=norm (v) ; 

%  angular  momentum 
h=cross (r , v) ; 

hmag=norm(h) ; 

%  eccentricity  vector 
e=cross (v, h) -r/rmag; 
ecc=norm{e) ; 

if  ecc<le-10,  %  Check  if  eccentricity  is  zero 

e  =  [0  0  0]  ; 


%  km 
%  sec 
%  km/sec 
%  km"3/sec^2 
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ecc=0 ; 

end 

%  node  vector 
n  =  cross (K,h) ; 

nmag=norm(n) ; 

if  nmag<le-10,  %  Check  if  inclination  is  zero 

n  =  [0  0  0]; 
nmag=0; 

end 

%  Find  the  true  anomaly  at  epoch 
nu  =  truanom(r , V, e,n) ; 

%  Compute  the  inclination 
inc  =  incl(h); 

%  Compute  the  longitude  of  ascending  node 
Omega=lan(n) ; 

%  Compute  the  argument  of  perigee 
omega  =  aop(n,e,r); 

%  Compute  the  rotation  matrix 
dcm=rot (Omega, inc, omega) ; 

direct=linspace (pi/2, 0,5) ; 
x=cos (direct) ;  y=sin (direct) ; 

2=2eros (1,5) ; 

for  k  =  1:5, 

dv=[x(k)  y(k)  2(k)]; 

if  abs  (x(k)  )  <le-‘12 
x(k) ^0; 

end 

if  abs (y (k) ) <16-12 
y(k)=0; 

end 

if  abs (y (k) ) <le-12 
y(k)=0; 

end 

end 

2=2eros (1,5) ; 
epoch=0 ; 

dvl  =  (dv_tot* [x(l)  y(l)  z(l)]'); 
dv2  =  (dv_tot* [x (2)  y(2)  z(2)]'); 
dv3  =  (dv_tot* [x(3)  y(3)  z(3)]'); 
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dv4  =  (dv_tot* [x (4 )  y{4)  z{4)]*); 
dv5  =  (dv_tot* [x (5)  y(5)  z(5)]’); 

rot  =  cz  (nu__new*dtr)  ; 

vnew_l  =  (rot* (vl  +  dvl 

vnew_2  =  (rot* (vl  +  dv2 

vnew_3  =  (rot* (vl  +  dv3 

vnew_4  =  (rot* (vl  +  dv4 

vnew_5  =  (rot* (vl  +  dv5 

%  Print  the  old  position  and  velocity  vectors  and  the  new  velocity 
vectors 

r_prev  =  rl; 
v_prev  =  vl; 

dv  =  [v_j)rev;  vnew_l ;  vnew_2 ;  vnew__3 ;  vnew_4 ;  vnew__5  ]  ; 

%fid=fopen( ' new_orb. vec ’ ,’w'); 

%fprintf (f id, ' ObjVectors  x  y  z  vx  vy 

vz  epoch\n'); 

%fprintf (fid, * %6 . Of %11 . 5f %11 . 5f%ll . 5f %10 . 6f%10 . 6f %10 . 6f %11 . 4f\n' , [0  r_prev 
v_prev  epoch] ) ; 

%fprintf (fid, ’ %6 . Of %11 . 5f%ll . 5f %11 . 5f %10 . 6f %10 . 6f %10 . 6f%ll . 4f \n ' , [1  r_prev 
vnew_l  epoch] ) ; 

%fprintf (fid, * %6 . Of %11 . 5f %11 . 5f %11 , 5f %10 , 6f%10 . 6f %10 . 6f %11 . 4f\n' , [2  r_prev 
vnew__2  epoch]  )  ; 

%fprintf (fid, * %6 . Of %11 . 5f %11 . 5f %11 . 5f %10 . 6f%10 . 6f %10 . 6f %11 . 4f\n ’ , [3  r_prev 
vnew_3  epoch] ) ; 

%f print f (fid, ' %6 . Of %11 . 5f %11 . 5f %11 . 5f %10 . 6f %10 . 6f %10 . 6f %11 . 4f \n ’  ,  [4  r_prev 
vnew_4  epoch] ) ; 

%fprintf (fid, ’ %6 . Of %11 . 5f%ll , 5f %11 . 5f %10 . 6f%10 . 6f%10 . 6f %11 . 4f\n’ , [5  r_prev 
vnew_5  epoch] ) ; 

%fprintf  ( ’  \nObjVectors  x  y  z  vx  vy 

vz  \n ' ) ; 

%fprintf ( ’ %6 . Of %11 . 5f %11 . 5f %11 . 5f %10 . 6f %10 . 6f %10 . 6f\n* , [0  r_prev  v_jDrev  ] ) ; 
%fprintf ( ’ %6 . Of %11 . 5f %11 . 5f %11 . 5f %10 . 6f %10 . 6f %10 . 6f\n* , [1  r_prev  vnew_l  ]); 
%fprintf  ( '  %6  .  Of  %11 . 5f  %11 . 5f  %11 . 5f  %10 . 6f  %10 . 6f  %10 .6f\n'  ,  [2  r_prev  vnew__2  ]  )  ; 
%fprintf  ( ’  %6  .  Of  %11 . 5f  %11 . 5f  %11 . 5f  %10 . 6f  %10 . 6f %10 . 6f\n*  ,  [3  r_prev  vnew__3  ]  )  ; 
%fprintf (’%6.0f%11.5f%11.5f%11.5f%10.6f%10.6f%10. 6f\n* ,  [4  r_prev  vnew_4  ] )  ; 
%fprintf ( ' %6 . Of %11 .5f %11 . 5f %11 . 5f %10 . 6f %10 . 6f %10 .6f\n* , [5  r_j)rev  vnew_5  ]); 

[el_prev,el_next,rv_next]  ^elements  (rl , dv, epoch, nu__new)  ; 

C.  DETERMINATION  OF  ORBITAL  ELEMENTS 

%  First  bring  in  the  r  and  v  vectors  that  you  want  to  convert . 

%  r  is  the  position  vector 

%  dv  is  a  5x3  matrix  with  the  velocity  vectors  for  the  various  delta  v’s 
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function  [el_prev, el_next , rv_next]  =  elements (r, dv, epoch, nu_new) 


rl=  r; 

Vl  =  dV; 

%  Constants 
rtd=180/pi; 
dtr=l/rtd; 

I  =  [1  0  0]  ; 

J  =  [0  10]; 

K  =  [0  0  1]; 

AU  =  1.4959965e8; 
TU  =  5.0226757e6; 
SU  =  29.784852; 
mu  =  1; 


%  Orbit  values : 

[m  n] =size (dv) ; 
for  s  =  l:m, 


%  Convert:  to  canonical  units 
r-ri*l86/AU; 
v=vl (s, : ) /SU; 

rmag=norm{r) ; 
vmag=norm(v) ; 

%  angular  momentum 
h=cross (r, v) ; 
hmag=norm(h) ; 

%  eccentricity  vector 
e=cross (v, h) -r/rmag; 
ecc_n=norm(e) ; 

if  ecc_n<le-10,  %  Check  if  eccentricity  is  zero 

e  =  [0  0  0]; 
ecc_n=0; 

end 

%  semi -major  axis 
p  =  hmag^2/mu; 
a  =  p/ (l-ecc_n^2) ; 

■'o  node  vector 
n  =  crosii?, (K,h)  ; 

nriiag=norm  (n)  ; 

if  nmc.g<le-10,  %  Check  if  inclination  is  zero 

n  =  [C  ■)  0]  ; 


%  km 
%  sec 
%  km/sec 
%  km^3/sec'^2 
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nmag=0 ; 
end 

%  Find  the  true  anomaly  at  epoch 
nu  =  truanom(r,v,e,n) ; 


%  Compute  the  inclination 
inc  =  incl (h) ; 

%  Compute  the  longitude  of  ascending  node 
Omega=lan (n) ; 

%  Compute  the  argument  of  perigee 
omega  =  aop(n,e,r); 

%  Compute  the  rotation  matrix 
dcm=rot (Omega, inc, omega) ; 

%  Compute  the  eccentric  anomaly  and  mean  anomaly 
E  =  acos (  (ecc_n  +  cos (nu*dtr) )  /  (1+  ecc_n*cos (nu*dtr) )  ); 

if  nu>180 

E=2*pi-E; 

end 

M  =  (E  -  ecc_n*sin (E)  )  ; 

H=:0; 

mu_s=l; 

n  =  sgrt  (mu_s/a^3)  ; 
global  t 

M_next  =  M  -  n*t; 

M_nextl  =  M_next  -  f ix (M_next/2/pi) *2*pi ; 
global  M_next  ecc_n 

%  Next  eccentric  anomaly 
E_next  =  fzero { ' funl ' ,M_next) ; 

%  Compute  next  true  anomaly 

nu__next  =  acos ( (ecc_n-cos (E_next) ) / (ecc_n*cos (E_next) -1) ) *rtd; 

if  sin (E_next ) <0 , 

nu_next =3  6  0 - nu_next ; 

end 

%  Compute  rmag  for  next  time 
rmag_next  =  a*  (l-ecc_n*cos  (E__next)  )  ; 

%  Compute  the  r/v  at  the  next  time 

[r_next  v_next  ]  =el ip2  ( rmag_next ,  nu_new ,  nu__next ,  ecc_n ,  a)  ; 
r__next=r_next  ; 
v__next  =  v__next  ; 

el_prev(s,  :)  =  [a  ecc_n  Omega  inc  omega  M  n  E]  ; 
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oV>  o\o  o\o 


el^next  (s ,  : )  =  [M_nextl  E_next  nu_next  rmag_next]  ; 
rv_next (s , : ) = [r_next  v_next] ; 

end 


D,  ELIP  1 


Calculation  of  the  Eliptical/circular  orbit 


LCDR  Wade  Knudson 

Naval  Postgraduate  School 

Orbit  Visualization  Routines 

Inputs:  rmag,dcm,nu 

User  input:  a,ecc,newinc 

Output:  r,v,  and  nu  for  the  new  orbit 


Files  called:  cx,cz 


[r,v,nu_new]  =  elip (rmag, dcm, nu) 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

★ 

* 

* 


function  [r, v, nu_new,dcm2 , ecc]  =  elipl (rmag, dcm, nu) 
%  Constants 


rtd 

-  180/pi; 

dtr 

=  1/rtd; 

AU  = 

1.4959965e8; 

% 

km 

TU  = 

5.0226757e6; 

% 

sec 

SU  = 

29.784852; 

“o 

km/sec 

mu  = 

1; 

% 

km^3/sec^2 

These  are  the  controlling  parameters  of  the  new  ellipse 


3'6'5'S'o'o'6'o^'S'5'5‘S'5'o'o‘b'o'b'o'6'5'S'S'€^'^'^'o'o*^'o'o'o'o<'5^'o'o<‘o'o'o'S^'o'o*^'o'o'o'o%^%^%^%^ 


a=l.l; 
ecc=0 . 1; 
r_p=a* (1-ecc) ; 

p=  r_jD*(l-ecc);  %AU  /  (1-ecc) /AU; 
newinc=0; 


if  ecc  ==  0, 

nu_new  =  0 ; 
a  =  rmag; 


%  If  circular  then  no 
%  rotation  to  argument 
%  of  periapsis  (nu__new) 
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p  =  a; 

else 

p  =a*  (l-ecc'^2)  ; 
arg= (p/rmag-1) /ecc; 

%if  arg-l<le-10 
%arg=l; 

%end 

nu_new=acos (arg) *rtd;  %  true  anomaly 

end 

rp  =  p/ (1+ecc) ; 
ra  =  p/ (l-ecc) ; 

%Perifocal  coordinates 
rl  =  trmag*cos (nu_new*dtr)  rmag*sin (nu_new*dtr)  0] ; 
vl  =  sqrt (mu/p) * [-sin (nu_new*dtr)  ecc+cos (nu_new*dtr)  0]  ; 

h=:cross  (rl, vl)  ; 
e=cross  (vl, h)  /mu-r/rmag; 


%  periapsis 
%  apoapsis 


%  rotate  the  new  orbit  to  the  line  of  node  (-nu) ,  then  to  desired 
inclination, 

%  and  finally  rotate  nu_2  to  get  epoch  to  r. 


dcm2=  dcm*cz ( -nu^dtr) *cx ( -newinc*dtr) *cz (nu_new*dtr) ; 
xyz_r=  (dcm2*rl '  )  '  ; 

r  =  xyz_r*AU/le6 ; 

V  =  (dcm2*vl' ) ' *SU; 


E.  ELIP  2 

% - 

% 

%  Calculation  of  the  Eliptical/circular  orbit 

% 

% 

%  LCDR  Wade  Knudson 

%  Naval  Postgraduate  School 

%  Orbit  Visualization  Routines 

% 

%  Input  s :  rmag , dcm , nu 

% 

%  Output:  r,v,  and  nu  for  the  new  orbit 
% 

%  Files  called:  cx,cz 

% 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


* 

* 
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[r,v,nu_new]  =  elip (rmag,dcm,nu) 


function  [r,v]  =  elip2 (rmag, dcm, nu_prev, ecc, a) 

%  Constants 
rtd=180/pi; 
dtr=l/rtd; 

I  -  [10  0]; 

J  =  [0  10]; 

K  =  [0  0  1]; 

AU  =  1.4959965e8;  %  km 

TU  =  5.0226757e6;  %  sec 

SU  =  29.784852;  %  km/sec 

mu  =  1;  %  km^3/sec^2 

p  =  a* {l-ecc^2)  ; 

=  P/ (l+ecc) ; 
ra  =  p/ (1-ecc) ; 

temp  -  cz  (nu__prev*dtr)  ; 

r  =  (temp* [rmag*cos (nu j)rev*dtr) ;  rmag*sin (nu_prev*dtr) ;  0] )  ' ; 

V  =  (temp* (sqrt (mu/p) ) * [-sin (nu_prev*dtr) ;  ecc+cos (nu_prev*dtr) ;  0] ) ’ 

r=r*AU/le6; 

v=v*SU; 


%  periapsis 
%  apoapsis 


F.  FUNCTION 

%  Solve  Kepler's  equation  for  E 
function  x  =  fun(E) 
global  M__prev  ecc__a 

X  =  Mj)rev  -  E  +  ecc_a*sin  (E)  ; 

G.  FUNCTION  1 

%  Solve  Kepler's  equation  for  E 
function  [x]  =  funl(E) 
global  M_next  ecc_n 


X  =  M_next  -  E  +  ecc_n*sin (E) ; 

aPARAM 

%  Solves  for  a, ecc, Omega, i, omega, nu,E,M,dcm,  and  rmag  from  r,v 
function  [a, ecc, O, inc, o, nu, E,M, dcm, rmag]  =  param(rl,vl) 

%  Constants 

rtd=lSC/pi; 

dtr=l/rtd; 
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I  =  [10  0]; 

J  =  [0  10]; 

K  =  [0  0  1]; 

AU  =  1.4959965e8;  %  km 
TU  =  5.0226757e6;  %  sec 
SU  =  29.784852;  %  km/sec 

mu  =  1;  %  km'^3/sec'*'2 

%  Earth  values : 

r=rl*le6/AU; 

v=vl/SU; 

rmag=norm(r) ; 
vmag=norm(v) ; 

%  angular  momentum 
h=cross (r, v) ; 

hmag=norm(h) ; 

%  eccentricity  vector 
e=cross (v, h) -r/rmag; 
ecc=norm(e) ; 

if  ecc<le-10,  %  Check  if  eccentricity  is  zero 

e  =  [0  0  0]; 
ecc=0; 

end 

%  node  vector 
n  =  cross (K,h) ; 

nmag=norm(n) ; 

if  nmag<le-10,  %  Check  if  inclination  is  zero 
n  -  [0  0  0]; 
nmag=0 ; 

end 

%  Find  the  true  anomaly  at  epoch 
nu  =  truanom (r , V, e, n) ; 

%  Compute  the  inclination 
inc  =  incl (h) ; 


p  =  hmag'^2/mu; 
a  =  p/ {l“ecc^2) ; 

%  Compute  the  longitude  of  ascending  node 
0=lan (n) ; 

%  Compute  the  argument  of  perigee 
o  =  aop (n, e,r) ; 

%  Compute  the  rotation  matrix 
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dcm=rot ( 0 , inc , o ) ; 


E  =  acos (  (ecc  +  cos(nu*dtr)  )  /  (1  +  ecc*cos {nu*dtr) )  ) ; 
M  =  E  -  ecc*sin(E); 
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720  Irwin  Ave.  (Stop  82-02) 

Falcon  Air  Force  Base 

Colorado  Springs,  Colorado  80912-7210 
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9.  Space  Warfare  Center/AE  1 

ATTN:  Deputy  Director 

720  Irwin  Ave.  (Stop  82-02) 

Falcon  Air  Force  Base 

Colorado  Springs,  Colorado  80912-7210 

10.  Space  Warfare  Center/XR  1 

ATTN;  Deputy  Director 

720  Irwin  Ave.  (Stop  82-02) 

Falcon  Air  Force  Base 

Colorado  Springs,  Colorado  80912-7210 

1 1 .  Lieutenant  Commander  Wade  E.  Knudson  2 
Weapons  Test  Squadron 

Pt.  Mugu,  California  93042-5001 
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